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ABSTRACT 

Stellar physics and evolution calculations enable a broad range of research 
in astrophysics. Modules for Experiments in Stellar Astrophysics (MESA) is a 
suite of open source, robust, efficient, thread-safe libraries for a wide range of ap- 
plications in computational stellar astrophysics. A 1-D stellar evolution module, 
MESA star, combines many of the numerical and physics modules for simulations 
of a wide range of stellar evolution scenarios ranging from very-low mass to mas- 
sive stars, including advanced evolutionary phases. MESA star solves the fully 
coupled structure and composition equations simultaneously. It uses adaptive 
mesh refinement and sophisticated timestep controls, and supports shared mem- 
ory parallelism based on OpenMP. State-of-the-art modules provide equation of 
state, opacity, nuclear reaction rates, element diffusion data, and atmosphere 
boundary conditions. Each module is constructed as a separate Fortran 95 li- 
brary with its own explicitly defined public interface to facilitate independent 
development. Several detailed examples indicate the extensive verification and 
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testing that is continuously performed, and demonstrate the wide range of capa- 
bilities that MESA possesses. These examples include evolutionary tracks of very 
low mass stars, brown dwarfs, and gas giant planets to very old ages; the complete 
evolutionary track of a 1M star from the pre-main sequence to a cooling white 
dwarf; the Solar sound speed profile; the evolution of intermediate mass stars 
through the He-core burning phase and thermal pulses on the He-shell burning 
AGB phase; the interior structure of slowly pulsating B Stars and Beta Cepheids; 
the complete evolutionary tracks of massive stars from the pre-main sequence to 
the onset of core collapse; mass transfer from stars undergoing Roche lobe over- 
flow; and the evolution of helium accretion onto a neutron star. MESA can be 
downloaded from the project web siteQ 

Subject headings: stars: general — stars: evolution — methods: numerical 
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1. Introduction 



Much of the information that astronomers use to study the universe comes from starlight. 
Interpretation of that starlight requires a detailed understanding of stellar astrophysics, 
especially as it relates to stellar atmospheres, structure, and evolution. Stellar structure and 
evolution models underpin much of modern astrophysics as they are used to analyze: the Sun 
through helioseismology (e.g., Bahcall et al.||1998 ), the pulsational properties of many nearby 
stars with asteroseismic data from, e.g., Corot ( Degroote et al.|2009 ) and Kepler (Gilliland et 
aLl[20T0| , the color -magnitude diagrams of resolved stellar and sub-stellar populations in the 



Milky Way and nearby galaxies (e.g., VandenBerg 2000 Dotter et al. 2010), the integrated 



light of distant galaxies and star clusters via population synthesis techniques (e.g., Worthey 
1994 Coelho et al. 2007), stellar yields and galactic chemical evolution (e.g., Timmes et 



al. 1995), physics of the first stars (Fujimoto et al. 2000), and a variety of aspects in time 
domain astrophysics (e.g., LSSlQ. 



Stellar evolution is broadly recognized as the first key problem in computational as- 
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trophysics. The introduction of electronic computers enabled the solution of the highly 
non-linear, coupled differential equations of stellar structure and evolution, and the first 



detailed reports of computer programs for stellar evolution soon appeared ( 


Iben & Ehrman 


1962 


Henyey et al. 


1964; 


Hofmeister et al. 


1964 


Kippenhahn et al. 


1967 


). Implicit in 



the development of these codes was a sufficiently mature theoretical understanding of stars 
( Chandrasekhar| 1938 Schwarzschild 1958, see as well the compilation of references later 
in this section), development of a concise yet sufficiently accurate treatment of convection 



( Bohm- Vitense 1958), as well as a better understanding of the properties of stellar matter, 
including nucleosynthesis ( Burbridge et al7||1957 Cameron|1957 ). Further improvements and 
alternative implementations became available addressing, for example, the numerical stabil- 



ity of computations (Sugimoto 1970), more efficient methods for following shell burning in 
low mass stars (Eg gleton|1971 ), and the hydrodynamics of advanced burning in massive stars 
(Weaver et al. 1978[ ). Progress continues on stellar evolution codes, with code developments 
and comparisons often facilitated by the opening of new observational windows. For example, 



the participants ( Christensen-Dalsgaard 


2008 DegPInnocenti et al. 


2008 Demarque et al. 


2008 Eggenberger et al. 2008 Hui-Bon-Hoa 


2008 


Morel & Lebreton 


2008; Roxburgh 2008 


Scuflaire et al. 2008; Ventura et al. 2008 


Weiss & Schlattl 2008) in 


the CoRoT Evolution 
/e sample of the active 


and Seismic Tools Activity (Lebreton et al. 


2008) 


are a represent at i 1 



community. 

Modules for Experiments in Stellar Astrophysics (MESA) began as an effort to improve 



upon the EZ stellar evolution code (Eggleton 1971 Paxton 2004). It employs modern soft 



ware engineering tools and techniques to target modern computer architectures that are 
significantly different from those available to the pioneers half a century ago. As the pieces 
of the new system started to emerge, it became clear that the parts would be of greater 
value than the whole if they were carefully structured for independent use. MESA includes a 
new 1-D stellar evolution code, MESA star, but is designed to be useful for a wide range of 
stellar physics applications. The physical inputs to stellar evolution models, like the equation 
of state, opacities, and nuclear reaction networks, have a broader application than stellar 
evolution calculations alone. MESA is designed so that each of the individual components is 
usable on its own, with the intention of facilitating verification test suites amongst different 
codes and encouraging new computational experiments in stellar astrophysics. 

MESA star approaches stellar physics, structure, and evolution with modern, sophisti- 
cated numerical methods and updated physics that give it a very wide range of applicability. 
The numerical and computational methods employed allow MESA star to consistently evolve 
stellar models through challenging phases, e.g., the He core flash in low mass stars and ad- 
vanced nuclear burning in massive stars, that have posed substantial challenges for stellar 
evolution codes in the past. 
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MESA is open source: anyone can download the source code, compile it, and run it for 
their own research or education purposes. It is meant to engage the broader community of 
astrophysicists in related fields and encourage contributions in the form of testing, finding 
and fixing bugs, adding new capabilities, and, generally, sharing experience with the MESA 
community. The philosophy and guidelines of MESA are described in more detail in the MESA 
manifesto (see Appendix |A|). 



This paper serves as an introduction to MESA and demonstrates its current capabili- 
ties. We assume that the reader is familiar with the basic stellar physics and numerical 
methods, both of which are essential to arrive at meaningful solutions when using MESA. 



For background material we refer the reader to Eddington 


(1926 


), Chandrasekhar (1939), 


Schwarzschild 


1958) 


, Cox & Giuli (1968), Clayton 


(1984 


), 1 


ben ( 


1991), Hansen & Kawaler 


(1995), Arnett 


(1996 


, and Kippenhahn & Weigert 


(1996 


• 



The MESA codebase is in constant development, and future capabilities and applications 
will be detailed in subsequent papers. The paper is outlined as follows: £j2] explains the 
design and implementation of MESA modules; QS] describe the numerical, microphysics, and 
macrophysical modules; ^describes the stellar evolution module MESA star; §[7] presents a 
series of tests and code comparisons that serve as rudimentary verification and demonstrates 
the broader capabilities of of MESA star; and $8] summarizes the material presented. 
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Table 1. Variable Index 



Name 



Description 



First appears 



A 
a 
a 

"MLT 
C 

C P 
C v 

Cs 
Xp 
AT 
D 

Dov 

A 

S 

dm 

dP 3 

dT s 

St 

E 



F 
f 
9 

r 

r :i 

v ad 

v rad 

V T 

Ks 

L 

Lconv 

A 

Xp 

m 

M c 
M m 

H 

Me 

N 

V 
P 



atomic mass number 
acceleration at the cell face 
order of convergence 
mixing length parameter 

"spacetime" parameter for convergence study 

specific heat at constant pressure 

specific heat at constant volume 

sound speed 

= dlnP/(ilnp| T 

= dlnP/dlnT| p 

Eulerian diffusion coefficient 

overshoot diffusion coefficient 

grid difference 

time difference 

mass of a cell 

P difference between surface and center of first cell 
T difference in between surface and center of first cell 
timestep 
energy 

nuclear energy generation in ergs/g 

power per unit mass (nuclear, thermal neutrino, gravity) 

Fermi energy 

mass flow rate 

overshoot mixing parameter 

local gravity 

= d\nP/d\np\ s 

Coulomb coupling parameter 

= dlnT/d\np\s + 1 

adiabatic temperature gradient 

radiative temperature gradient 

actual temperature gradient 

opacity at the surface of the outermost cell 

total luminosity 

convective luminosity 

mixing length (a MLT A P ) 

pressure scale height 

mass interior to cell 

inner mass (not modeled) for central BC 
modeled mass 

mean molecular weight per gas particle 
mean molecular weight per electron 
Brunt-Vaisala frequency 

dimensionless electron degeneracy parameter 
total pressure 
gas pressure 

pressure at surface of outermost cell 



(4.4 



(6.2 



E 6.7 



in 



(6.7 



(1.2 



(4.2 



(7.2.2 



(4.2 



(4.2 



(5.2 



(5.2 



(6.5 



(6.5 



(6.2 



(6.2 



(6.2 



(6.4 



(4.2 



(4.5 



(6.2 



( 7.1 



(6.2 



(5.2 



(5.1 



(4.2 



(1.2 



(4.2 



(4.2 



(5.1 



(5.1 



(5.3 



(5.1 



(5.1 



(5.1 



(5.1 



(6.2 



(6.6 



(6.6 



(4.2 



(1.2 



(7.2.2 



(6.5 



(4.2 



(4.2 



( 5.3 
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Table 1 — Continued 



Name 



Description 



First appears 



q relative mass coordinate 

p density 

R total radius 

-Rcz radius of the base of the solar convective zone 

r radius at the cell face 

S entropy 

a Lagrangian diffusion coefficient 

T temperature 

T s temperature at surface of outermost cell 

T c ff effective temperature 

t s optical depth at the surface of the outermost cell 

t optical depth 

v velocity at the cell face 

v c timestep control target 

fconv convective velocity 

v t timestep control variable 

w diffusion velocity 

X H mass fraction 

Xi mass fraction of the i th isotope 

^ relative difference in convergence study 

Y He mass fraction 

Y e electrons per baryon (Z/A) 

Z metals mass fraction (1 — X — V) 

Z atomic number 

2 distance from convective boundary 



[6^ 
j O 
[53 

!5T 
j O 
iO 
[53 
[O 
iO 
\ K2 

j O 
j 5T 

; 573 

j 672" 
[67? 
[ O 
[773 

[ 472 

[ 574 

f B~2 
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2. Module design and implementation 

Each MESA module is responsible for a different aspect of numerics or physics required 
to construct computational models for stellar astrophysics. Each has a similar organization: 
a public interface, a private implementation, a makefile to create a library, and a test suite 
for verification. Each module includes an installation script that builds the library, tests it, 
and, if the test succeeds, exports it to the MESA libraries directory. Comparisons between 
local and expected results are carried out with the open source ndif f utility]^] There is a 
global install script for MESA that performs the installation of each of the modules in the 
required order to satisfy dependencies. The installation on UNIX-like systems, including 
Linux and Mac OS X requires a modern, up-to-date Fortran compiler^] A template module, 
package_template, exists for initiation of new modules by the community. All current MESA 
modules are listed in Table [2j along with the function they perform and the section in this 
paper where the description resides. 

The MESA modules are "thread-safe" — meaning that more than one process can execute 
the module routines at the same time — allowing applications to utilize multicore processors. 
A module is thread-safe if all of its shared data is read-only after initialization. This prohibits 
the use of common blocks and "SAVE" statements. Working memory must be allocated 
as local variables of routines or allocated dynamically. To take full advantage of shared 
memory on multicores, an operation that is performed in parallel needs to fit in the processor 
cache. Evaluations of local microphysics, such as the equation of state, opacity, and nuclear 
reaction networks can be carried out in parallel using the OpenMP application programming 
interface]^] The capability of MESA star to take advantage of multithreading is discussed in 

go 



3. Numerical methods 

MESA includes several modules that provide numerical methods. The following briefly 
describes each one presently available and references the relevant literature (or web-based 
resource) where the full description resides. 



5 See 



http : //www . math . Utah . edu/~beebe/ sof tware/ndif f / 



MESA installs its own copy of ndiff the 



first time the main installation process is performed. 

information about supported compilers and installation is provided on the MESA project website. 



E http : / /openmp . org 
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Table 2. MESA Module Definitions and Purposes 



Name 



Type 



Purpose 



Section 



alert 


utility 


error handling 


atm 


microphysics 


grey and non-grey atmospheres; tables and integration 


const 


utility 


numerical and physical constants 


chem 


microphysics 


properties of elements and isotopes 


diffusion 


macrophysics 


gravitational settling and chemical and thermal diffusion 


eos 


microphysics 


equation of state 


interp_ld 


numerics 


1-D interpolation routines 


interp_2d 


numerics 


2-D interpolation routines 


ionization 


microphysics 


average ionic charges for diffusion 


jina 


macrophysics 


large nuclear reaction nets using reaclib 


kap 


microphysics 


opacities 


karo 


microphysics 


alternative low-T opacities for C and N enhanced material 


mlt 


macrophysics 


mixing length theory 


mtx 


numerics 


linear algebra matrix solvers 


net 


macrophysics 


small nuclear reaction nets optimized for performance 


neu 


microphysics 


thermal neutrino rates 


num 


numerics 


solvers for ordinary differential and differential-algebraic equations 


package_template 


utility 


template for creating a new MESA module 


rates 


microphysics 


nuclear reaction rates 


screen 


microphysics 


nuclear reaction screening 


star 


evolution 


1-D stellar evolution 


utils 


utility 


miscellaneous utilities 


weaklib 


microphysics 


rates for weak nuclear reactions 



3 

^3 
IT 

TJ 
El 

4-1 

H 
= 

57T 

I 

C5 

4~5 

TIT 

2 
4.4 
43 

Wi 

3 

4.5 
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The mtx module provides an interface to linear algebra routines for matrix manipulation. 
A large set of BLAS and LAPACK routines are included, but the mtx module can easily be 
modified to accept these routines from other linear algebra packages, e.g. GotoBLAS0 or 
the Intel Math Kernel Librarjj^] (MKL). Sparse matrix operations are supported, including a 
subset of the SPARSKIT sparse matrix iterative solveij^] and an interface to the Intel version 
of the PARDISO sparse matrix direct solver. The routines in num make use of these matrix 
routines. 

Modules interp_ld and interp_2d deal with 1-D and 2-D interpolation, respectively. 
One dimensional interpolation is carried out using either a piecewise monotonic cubic method 



(Huynh 1993 Suresh & Huynh 1997) or a monotonicity-preserving method (Steffen 1990). 



Compared to the piecewise monotonic method, the monotonicity-preserving method is stricter 



and does not allow an interpolated value to range outside of the given values (Steffen 1990). 



Module interp_2d includes parts of the PSPLINE package^] and routines by both Akima 



(1996) and Renka (1999) for bivariate interpolation and surface fitting on a grid or with a 



scattered set of data points. Both single- and double-precision versions of the 2-D interpo- 
lation routines are provided. 

Module num provides a variety of solvers for stiff and non-stiff systems of ordinary 
differential equations (ODEs) and a Newton-Raphson solver for multidimensional, nonlinear 
root-finding. The family of ODE solvers is derived from the routines of |Hairer fc Wanner 



(1996). The non-stiff ODE class are explicit Runge-Kutta integrators of orders 5 and 8 



with dense output, automatic stepsize control, and optional monitoring for stiffness. The 
stiff ODE solvers are linearly implicit Runge-Kutta, with 2nd, 3rd, and 4th order versions 
and two implicit extrapolation integrators of variable order: either midpoint or Euler. All 
integrators support dense, banded, or sparse matrix routines, analytic or numerical difference 
Jacobians, explicit or implicit ODE systems, dense output, and automatic stepsize control. 

The Newton-Raphson solver for multidimensional, nonlinear root-finding supports square, 
banded, and sparse matrices and analytic or automatic numerical differencing for the Ja- 
cobians. It has the ability to reuse Jacobians and employs a line search method to give 
improved convergence. The multidimensional Newton-Raphson solver is used by MESA star 
to solve highly non-linear systems of differential-algebraic equations with tens of thousands 



e http : //www . tacc . utexas . edu/ resources/ software 




'http : / / software . intel . com/ en-us/ intel-mkl/ 






6 http : //www-users . cs.umn.edu/saad/ sof tware/SPARSKIT/ sparskit .html 


s http : //w3 . pppl . gov/NTCC/PSPLINE 
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of variables (see §|6j). The structure of the Newton- Raphson solver is derived from Lesaffre's 
version of the Eggleton stellar evolution code ( Eggleton||1971 Pols et aL]|1995 Lesaffre et al. 
2006) and some details of the implementation will be described in 



The alert module provides a framework for reporting messages, including errors, to 
the terminal. The utils module provides a number of functions for checking if a variable 
has been assigned a bad value (e.g., NaN or Infinity) and tracking Fortran I/O unit numbers 
in use. It also provides subroutines for basic file I/O and for allocating arrays of different 
types and dimensions, including a Fortran implementation of a hash tree that is used by the 
stellar evolution module to update the model mesh. Programs and scripts that are used for 
testing that each module has compiled correctly are stored in utils. 



4. Microphysics 

The MESA microphysics modules provide the physical properties of stellar matter, with 
each module focusing on a separate aspect of the physics. 



4.1. Mathematical constants, physical and astronomical data 

The MESA module const contains mathematical, physical, and astronomical constants 
relevant to stellar astrophysics in cgs units. The primary source for physical constants is the 
CODATA Recommended Values of the Fundamental Physical Constants ( Mohr et al.| [2008). 



Values for the Solar age, mass, radius, and luminosity are taken from Bahcall et al. (2005) 



The MESA module chem is a collection of data, functions, and subroutines to manage 
the chemical elements and their isotopes. It contains basic information about the chemi- 
cal elements and their isotopes from Hydrogen through Uranium. It includes routines for 
translating between atomic weights and numbers and isotope names. It contains full list- 
ings of Solar abundances on several scales ( Anders fc Greves"se||1989 Grevesse fc Noels||1993 



Grevesse fc Sauva"I||1998 Lodders||2003 Asplund et al.||2004 ). Module chem contains a frame- 
work for the user to provide an arbitrary set of species in a text file. 



4.2. Equation of state 



The equation of state (EOS) is delivered by the eos module. It works with density, p, and 
temperature, T, as independent variables. These are the natural variables in a Helmholtz 



13 



free energy formulation of the thermodynamics. However, as some calculations are more 
naturally performed using pressure, P, and T (as in a Gibbs free energy formulation), a 
simple root find can provide p given the desired P gas = P — aT 4 /3 and T. While conceptually 
simple, this can impose a substantial computational overhead if done for each eos call. To 
alleviate this computational burden, the root finds are pre-processed, creating a set of tables 
indexed by P gas and T. As a result, the runtime cost of evaluating eos using P gas and T is 
the same as for using p and T, as long as the P gas — T requests are within the pre-computed 
ranges. When outside those ranges, the root find is performed during runtime, slowing the 
computations. 



The MESA p-T tables are based on the 2005 update of the OPAL EOS tables (Rogers 



& Nayfonov 2002). To extend to lower temperatures and densities, we use the SCVH ta- 



bles (Saumon et al. 1995), and construct a smooth transition between these tables in the 
overlapping region that we define (shown by the blue dotted lines in Figure [T| . The lim- 
ited thermodynamic information available from these EOSs restricts our blending to the 
output quantities listed in Table |3j The resulting MESA tables are more finely gridded than 
the original tables (so that no information is lost) and are provided at six X and three Z 
values: X = (0.0,0.2,0.4,0.6,0.8,1.0) and Z = (0.0,0.02,0.04) in keeping with the OPAL 
tables, allowing for Helium rich compositions. In order to save space, the MESA tables are 
not rectangular in the independent variables. Instead, the region occupied by usual stel- 
lar models is roughly rectangular in the stellar modeling motivated variables, logT and 
\ogQ = logp — 21ogT + 12. The range in logT is from 2.1 to 8.2 in steps of 0.02 and the 
range in log Q is from -10.0 to 5.69 in steps of 0.03. Partials with logT and log Q are derived 
from the interpolating polynomials, while partials with respect to log p then follow. The 
resulting region of these MESA tables is that inside of the dashed black lines of Figure [TJ The 
MESA P gas — T tables are rectangular in log T and log W = log P gas — 4 log T over a range 
-17.2 < log W < -2.9, and 2.1 < logT < 8.2. 



Outside the region covered by the MESA tables, the HELM ( Timmes fc Swesty||2000 ) and 



PC ( |Potekhin fc ChabrierfMol EOSs are employed. Both HELM and PC assume complete 



ionization and were explicitly constructed from a free energy approach, guaranteeing thermo- 
dynamic consistency. In nearly all cases, the full ionization assumption is appropriate since 
the OPAL and SCVH tables are used at those cooler temperatures where partial ionization 



is significant. Since the MESA tables are only constructed for Z < 0.04, eos uses HELM 



and PC for Z > 0.04 in the whole p — T plane. 



10 



We discuss the ionization states of trace heavy elements in I 



5.4 



-14- 




Fig. 1. — The p — T coverage of the equations of state used by the eos module for Z < 0.04. 
Inside the region bounded by the black dashed lines we use MESA EOS tables that were 
constructed from the OPAL and SCVH tables. The OPAL and SCVH tables were blended 
in the region shown by the blue dotted lines, as described in the text. Regions outside of 
the black dashed lines utilize the HELM and PC EOSs, which, respectively, incorporate 
electron-positron pairs at high temperatures and crystallization at low temperatures. The 
blending of the MESA table and the HELM/PC results occurs between the black dashed lines 
and is described in the text. The dotted red line shows where the number of electrons per 
baryon has doubled due to pair production, and the region to the left of the dashed red line 
has Ti < 4/3. The very low density cold region in the leftmost part of the figure is treated 
as an ideal, neutral gas. The region below the black dashed line labeled as T = 175 would 
be in a crystalline state for a plasma of pure oxygen and is fully handled by the PC EOS. 
The red dot-dashed line shows where MESA blends the PC and HELM EOSs. The green 
lines show stellar profiles for a main sequence star (M = l.OM ), a contracting object of 
M = O.OO1M and a cooling white dwarf of M = O.8M . The heavy dark line is an evolved 
25M star that has a maximum infalling speed of 1000 km s _1 . The jagged behavior reflects 
the distinct burning shells. 
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Table 3. eos output quantities and units 



Output Definition Units 



p 

1 gas 


gas pressure 


ergs cm~ 3 




E 


internal energy 


ergs g -1 




S 


entropy per gram 


ergs g _1 PC 


-i 


dE/dp\ T 




ergs cm 3 g~ 


-2 


C v 


specific heat at constant V = l/p 


ergs g^ 1 K 


-1 


dS/dp\T 




ergs cm 3 g~ 2 


K 


dS/dT\ p 




ergs g^ 1 K 


-2 


X P 


= d\nP/d\np\ T 


none 




XT 


= d\nP/tMT\ p 


none 




c P 


specific heat at constant pressure 


ergs g _1 K 


-1 


v ad 


adiabatic T gradient with pressure 


none 






= d\nP/d\np\s 


none 




r 3 


= d\nT/dlnp\s + 1 


none 




V 


ratio of electron chemical potential to ksT 


none 




n 


mean molecular weight per gas particle 


none 






mean number of free electrons per nucleon 


none 
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HELM was constructed for high temperatures (up to logT = 13) and densities (up 
to logp = 15), and accounts for the onset of electron-positron pair production at high 
temperatures. The dotted red line in Figure [T] shows where the number of electrons per 
baryon has doubled due to pair production. The domination of pairs in the plasma creates 
a region where Ti < 4/3 (to the left of the dashed red line). The blending region to HELM 
(from the MESA tables) is shown by the black dashed lines in Figure [TJ and can by modified 
by the user. In this transition region, the blend of the two EOSs is performed in a way that 
preserves thermodynamic consistency. Therefore, if each separate EOS satisfies Maxwell's 
relations, the blend will also satisfy them. To accomplish this, we linearly sum the EOS 
quantities Qi (i.e. P, E, S and their partial derivatives with respect to p and T) needed to 



satisfy Maxwell's relations (Timmes & Swesty 2000hPJ The blend is calculated by defining 



the boundary limits, inside of which we define a fractional "distance" , F, from the boundary. 
As F varies from zero to one, we use the smoothing function S = (1 — cos(F7r))/2 and for 
each of the nine quantities we construct Qi = SQf + (1 — S)Qf , where Qf and Qf are the 
outputs from the two EOSs. We then use these to rederive the thermodynamic quantities 
{Xp, Xt, C P , V ad , r 3 , Ti) delivered by the eos routine. 

In late stages of the cooling of white dwarfs, the ions in the core will crystallize. For pure 
oxygen, the crystallization limit corresponds to a value of the Coulomb coupling parameter, 
T w 175, shown by the black dashed line in Figure [1} In this region, we use the PC EOS, 
which accounts for the modified thermodynamics of a crystal, as well as carefully handling 
mixtures (e.g. carbon and oxygen). The blend between PC and HELM (as shown by the 
dot-dashed red lines in Figure [T]) is performed in the same manner as described above. In 
the dense liquid realm, the blending region is defined by the Coulomb coupling parameter, 
r = Z e 2 /aikpT, where aj is the mean ion spacing, and Z is the average ion charge. The 
default choice is PC for T > 80 and HELM for T < 40. The PC EOS is not constructed 
for arbitrarily low densities, forcing a transition to HELM at logp < 2.8, with the blend 
beginning at logp = 3.7. These boundaries may be re-defined by the user if needed. 

In addition to the two independent variables, the eos module requires as input X, Z, 
A (the mass-averaged atomic weight of metals), and Z (the mass-averaged atomic charge 
of metals). When operating in the regime where the PC EOS is implemented, the mass 
fractions for all isotopes with mass fractions above a specified minimum are needed (default 
is 0.01), allowing PC to correctly handle isotope mixtures. It returns a total of sixteen quan- 
tities (listed in Table [3j) as well as the partial derivatives of each quantity with respect to the 
independent variables. The tables are interpolated in the independent variables using bicu- 



11 The more conventional forms of these nine thermodynamic quantities are displayed in the first nine rows 
in Table [3] 
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bic splines from interp_2d with partial derivatives determined from the splines. Separate 
quadratic interpolations are performed in X and Z. 

The construction of eos tables as outlined above is the default option for MESA but the 
eos module has the flexibility to accept tables from any source so long as the tables conform 
to the MESA standard format. For example, the comparison with the Stellar Code Calibration 
project (Weiss et al. 2007) described in £7.1.2 utilizes tables constructed using FreeEOSQ 
The FreeEOS code does not cover the same range of p and T as SCVH+OPAL+HELM but 
the eos module is designed with this flexibility in mind: the table dimensions are specified 
in the table headers and the module dynamically allocates arrays of the appropriate size to 
hold them when the tables are read in. 



Since not all EOS sources may be in the tabular form desired by eos, we have created 
a module, other_eos, that provides the user an opportunity to incorporate their own EOS 
and use it with the stellar evolution module MESA star. 



4.3. Opacities 



The pre-processor make_kap resides within the kap module and constructs the MESA 
opacity tables by combining radiative opacities with the electron conduction opacities from 



Cassisi et al. (2007). In the rare circumstances where p or T are outside the region covered 



by Cassisi et al. (2007) (-6 < logp < 9.75 and 3 < logT < 9), the Iben (1975) fit to 



the Hubbard & Lampe (1969) electron conduction opacity is used for non-degenerate cases 



while the Yakovlev & Urpin (1980) fits are used for degenerate cases. Radiative opacities 



are taken from Ferguson et al. (2005) for 2.7 < logT < 4.5 and OPAL (Iglesias & Rogers 



1993, 1996) for 3.75 < logT < 8.7. The low T opacities of Ferguson et al. (2005) include 



the effects of molecules and grains on the radiative opacity. Tables from OP (Seaton et al. 



2005) can be used in place of OPAL as the table format is identical. The radiative opacity 



is dominated by Compton scattering for logT > 8.7 and is calculated using the equations of 
Buchler & Yuehl (fl976l up to a density of 10 6 g cm" 3 . We use the HELM EOS to calculate 



the increasing number of electrons and positrons per baryon when pair production becomes 
prevalent, an important opacity enhancement. 



The OPAL tables with fixed metal distributions are called Type 1 (Iglesias & Rogers 



1993, 1996) and cover the region 0.0 < X < 1 - Z and 0.0 < Z < 0.1. Additionally, there is 



support for the OPAL Type 2 ( Iglesias fc Rogers]|1996 ) tables that allow for varying amounts 



http : //f reeeos . sourcef orge . net 
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of C and beyond that accounted for by Z\ these are needed during helium burning and 
beyond. These have a range 0.0 < X < 0.7, 0.0 < Z < 0.1. 

The resulting kap tables cover the large range 2.7 < logT < 10.3 and —8 < logi? < 8 
(i? = p/Tq, so logi? = logp — 31ogT + 18), as shown by the heavy orange and black lines 
in Figure [2j The MESA release includes MESA opacity tables derived from Type 1 and 2 



OPAL tables, tables from OP, and Ferguson et al. (2005). The heavy orange lines delineate 



the boundaries where we use existing tables to make the MESA opacity table. The blended 
regions in Figure [2] are where two distinct sources of radiative opacities exist for the same 
parameters, requiring a smoothing function that blends them in a manner adequate for 
derivatives. The blend is calculated at a fixed logi? by defining the upper (logTJ/) and 
lower (\ogT L ) boundaries of the blending region in logT space, where Ku (k l ) is the opacity 
source above (below) the blend. We perform the interpolation by defining F = (logT — 
log Tl) / (log Tjj — logTx), and using a smooth function S = (1 — cos(F7r))/2 for 

log k = S log Ku (R, T) + (1 - S) log k l (R, T). (1) 

At high temperatures, the blend from Compton to OPAL (or OP) has logTu = 8.7 and 



\ogTi = 8.2. At low temperatures, the blend between Ferguson et al. (2005) and OPAL has 
logTJy = 4.5 and \ogT L = 3.75. 

The absence of tabulated radiative opacities for logi? > 1 and logT < 8.2 (the region 
below the heavy dashed line in Figure [2j leads us to use the radiative opacity at log R = 1 
(for a specific logT) when combining with the electron conduction opacities. This introduces 
errors in the MESA opacity table between log R = 1 and the region to the right of the dashed 
blue line in Figure [2] where conductive opacities become dominant. However, as we show in 
Figure |3j main sequence stars are always efficiently convective in this region of parameter 
space, alleviating the issue. 

The module kap gives the user the resulting opacities by interpolating in logT and 
logi? with bicubic splines from interp_2d. The user has the option of either linear or cubic 
interpolation in X and Z and can specify whether to use the fixed metal (Type 1) tables or 
the varying C and O (Type 2) tables. In the latter case, the user must specify the reference 
C and O mass fractions, usually corresponding to the C and O in the initial composition. 
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Fig. 2. — The sources of the standard MESA opacity tables. Construction of opacity tables 
requires incorporating different sources, denoted by the labels. The heavy orange lines 
denote regions where input tables exist for radiative opacities, whereas the heavy black 
lines extend into regions where we use algorithms to derive the total opacities, described 
in the text. Above the dashed red line, the number of electrons and positrons from pair 
production exceeds the number of electrons from ionization, and is accounted for in the 
opacity table. The opacity in the region to the right of the dashed blue line is dominated 
by electron conduction. Also shown are stellar profiles for stars on main sequence (M = 
0.1, 1.0, & IOOMq) or just below (a contracting M = O.O1M brown dwarf). 
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For requests outside the logT and logi? boundaries, the following is done. The region 
to the left of logi? = —8 and below logT = 8.7 is electron scattering dominated, so the 
cross-section per electron is density independent. However, the increasing importance of 
the Compton effect as the temperature increases (which is incorporated in the OPAL/OP 
tabulated opacities) must be included, so we use the opacity from the table at logi? = —8 
at the appropriate value of logT. For higher temperatures (logT > 8.7) electron-positron 
pairs become prevalent, as exhibited by the red dashed line that shows where the number of 
positrons and electrons from pair production exceeds the number of electrons from ionization. 
MESA incorporates the enhancement to the opacity from these increasing numbers of leptons 
per baryon. 

At the end of a star's life, low enough entropies can be reached that an opacity for 
logi? > 8 is needed. When kap is called in this region, we simply use the value at logi? = 8 
for the same logT. For regions where Z > 0.1, the table at Z = 0.1 is used. 

The resulting opacities for Z = 0.019 and Y = 0.275 are shown in Figure |3j both 
as a color code, and as contours relative to the electron scattering opacity, k = 0.2(1 + 
X) cm 2 g _1 . The orange lines show (top to bottom) where logi? = —8, logi? = 1 and 
logi? = 8. We show a few stellar profiles for main sequence stars as marked. The green 
parts of the line are where heat transfer is dominated by heat transport, requiring an opacity, 
whereas the light blue parts of the line are where the model is convective. As is evident, 
nearly all of the stellar cases of interest (shown by the green-blue lines) are safely within the 
boundaries or the MESA tables. The lack of radiative opacities in the higher density region to 
the right of log i? = 1 implies opacity uncertainties until the dark blue line is reached (where 
the conductive opacity takes over). However, the stellar models are convectively efficient in 
this region, so that the poor value for k does not impact the result as long as the convective 
zone's existence is independent of the opacity (the typical case for these stars, where the 
ionization zone causes the convection). 
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Fig. 3. — The resulting MESA opacities for Z = 0.019, Y = 0.275. The underlying shades 
show the value of k, whereas the contours are in units of the electron scattering opacity, 
no = 0.2(1 + X) cm 2 g _1 . The orange lines show (top to bottom) where logi? = —8, 
logi? = 1 and logi? = 8. Stellar interior profiles for main sequence stars of mass M = 
0.1,0.3,1.0,3.0 & 100M Q are shown by the green(radiative regions )-light blue(convective 
regions) lines. Electron conduction dominates the opacity to the right of the dark blue line 
(which is where the radiative opacity equals the conductive opacity). 

It is also possible to generate a new set of kap readable opacity tables using the make_kap 
pre-processor. The requirements are high-temperature radiative opacities in the standard 
OPAL format and low-temperature radiative opacities in the number and format provided 
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by Ferguson et al. (120051) J 13 Specific high-temperature radiative opacities can be made by 



using the OPAL site or the Opacity Project site 



15 



Since not all opacity sources can be placed in the tabular form desired by kap, we have 
created a module, other_kap, that provides the user an opportunity to incorporate their own 
opacity source. A simple flag tells MESA star to call other_kap rather than kap, allowing 
for experiments with new opacity schemes and physics updates. The first example of such 
an implementation that has now become a MESA module is karo. It was developed to study 
the stellar evolution effects of dust-driven winds in Carbon-rich stars, using the Rosseland 



opacities of Lederer & Aringer (2009) and the hydro-dynamical wind models of Mattsson et 

aDpnol). 



4.4. Thermonuclear and weak reactions 



The rates module contains thermonuclear reaction rates from Caughlan & Fowler (1988 



CF88) and |Angulo et al.| p999| NACRE), with preference given to the NACRE rate when 
available. The reaction rate library includes more than 300 rates for elements up to Nickel, 
and includes the weak reactions needed for Hydrogen burning (e.g. positron emissions, 
electron captures) , as well as neutron-proton conversions and a few other electron and neu- 
tron capture reactions. Significant updates to the NACRE rates have been included for 



ii 



N(p,7) 15 ( |Imbriani et aL]|2004[ ), triple-a ( |Fynbo et aL]|2005[ ), 14 N(a,7) 18 F flGorres et al 
2000) and 12 C(a,7) 16 (Kunz et al. 2002). In these special cases, the rate can be selected 



from CF88, NACRE, or the newer reference by the user at run time. 

The weaklib module calculates lepton captures and /3-decay rates for the high densities 
and temperatures encountered in late stages of stellar evolution. The rates are based on 



the tabulations of Fuller et al. (1985), Oda et al. (1994), and Langanke & Martmez-Pinedo 



(2000) for isotopes with 45 < A < 65. The most recent tabulations of Langanke & Martmez- 



Pinedo (2000) take precedence, followed by Oda et al. (1994), then Fuller et al. (1985). The 



user can override this to create tables using any combination of these or other sources. 

The screen module calculates electron screening factors for thermonuclear reactions in 



both the weak and strong regimes. The treatment has two options. One is based on Dewitt 



13 


http : //webs . wichita. edu/physics/opacity 
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et al. (1973) and Graboske et al. (1973). The othei 16 combines Graboske et al. (1973) in 



the weak regime and Alastuey & Jancovici (1978) with plasma parameters from Itoh et al. 



(1979) in the strong regime. 



The neu module calculates energy loss rates and their derivatives from neutrinos gen- 
erated by a range of processes including plasmon decay, pair annihilation, Bremsstrahlung, 
recombination and photo-neutrinos (i.e. neutrino pair production in Compton scattering). It 



is based on the publicly available routine (see footnote 16 ) derived from the fitting formulas 



of Itoh et al. (1996) 



4.5. Nuclear reaction networks 

The net module implements nuclear reaction networks and is derived from publicly 
available code (see footnote 16). It includes a "basic" network of 8 isotopes: 1 H, 3 He, 
4 He, 12 C, 14 N, 16 0, 20 Ne, and 24 Mg, and extended networks for more detailed calculations 
including coverage of hot CNO reactions, a-capture chains, (a,p) + (p,7) reactions, and heavy- 



ion reactions (Timmes 1999). In addition to using existing networks, the user can create 



a new network by listing the desired isotopes and reactions in a data file that is read at 
run time. The amount of heat deposited in the plasma by reactions is derived from the 



nuclear masses in chem, taken from the JINA Reaclib database (Rauscher & Thielemann 



2000 Sakharuk et al. 2006 Cyburt et al. 2010), and accounts for positron annihilations 



and energy lost to weak neutrinos, using Bahcall (1997, 2002) for the hydrogen burning 



reactions. The list of approximately 350 reactions is stored in a data file that catalogs the 
reaction name, the input and output species, and their heat release. 

The jina module is an alternative nuclear network module that specializes in large 
networks. It is based on the 'netjina' package by Ed Brown and uses the JINA Reaclib 



it 
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Table 4. Comparison of 1-zone Solar burn results at 10 Gyr 



Network log 10 e nuc login X( 1 H) log 10 X( 4 He) log 10 X( 12 C) log 10 X( 14 N) log 10 X( 16 0) 



jina 25 18.63757961 -3.87550319 -0.008144854 -4.40235799 -1.9195882 -3.07400339 
net 25 18.63685339 -3.87550517 -0.008145036 -4.40235799 -1.9195882 -3.07400333 
net 8 18.63675658 -3.93650004 -0.008137607 -4.39650625 -1.9135911 -3.04585377 
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Fig. 4. — A 1-zone hydrogen burn at constant T = 19 x 10 6 K and p = 100 g cm -3 by three 
different networks. The number following net or jina indicates the number of isotopes 
considered in that network. The 25 isotope networks expand on the 8 isotope network by 
including minor contributors to the pp and CNO cycles. The plot shows the evolution of the 
mass fraction abundances of the 10 most abundant isotopes and net energy generation per 
unit mass, e nuc (ergs g -1 ), as a function of time. The left-hand axis shows the mass fraction 
while the right-hand axis shows the net energy generation per unit mass. 



Table 5. Comparison of 1-zone He-burn results at 10 Gyr 



Network log 10 e nuc log 10 X( 12 C) log 10 X( 16 0) log 10 X( 22 Ne) log 10 X( 26 Mg) 



jina 200 17.9085633 -0.721578469 -0.108630252 -1.50380756 -4.01520633 
net 11 17.9086380 -0.721576540 -0.108630957 -1.50385214 -3.99780015 
net 8 17.9083877 -0.718866029 -0.107692784 
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Fig. 5. — Equivalent to Figure [4] but now showing a 1-zone helium burn at constant logT = 
8.1, logp = 4.0. The "net 11" network adds 18 0, 22 Ne, and 26 Mg to the 8 isotope network; 
"jina 200" includes about 200 isotopes up to 71 Ge. 
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database for thermonuclear reaction rates (Rauscher & Thielemann 2000; Sakharuk et al. 
2006 jF 7 ] the rates and weaklib modules for weak interactions, and screen for electron 



screening. Most importantly, it allows the user to create large nuclear networks by specifying 
the list of isotopes to consider. All nuclear reactions (both strong and weak) linking the 
isotopes in the set are automatically included in the network. In all, j ina covers more than 
76,000 nuclear reactions involving more than 4,500 isotopes. The jina module is slower 
than net for small networks but the flexibility and capacity to handle large networks make 
it advantageous in some cases. 

Both net and j ina include one-zone burn routines that operate on a user-defined initial 
composition, nuclear network, and a trajectory comprising density and temperature as a 
function of time. The one-zone burn routines interface with mtx and num, enabling the 
use of the sparse matrix solver, which substantially improves performance compared to the 
dense matrix solver for networks of more than a few hundred isotopes. Figures [4] and [5] 
demonstrate these one- zone routines operating on conditions appropriate for the Sun on the 
main sequence and for a core He-burning star, respectively. Both examples were evolved at 
fixed density and temperature for 10 Gyr. Each figure compares three networks of varying 
size in terms of the mass fractions and net energy generation per unit mass, e nuc (ergs g _1 ), 
produced. Tables [4] and [5] complement Figures [4] and |5j respectively, by listing the final values 
from each of the one-zone burn simulations. These comparisons indicate that the 8 isotope 
network produces results that agree with larger networks to 4-5 significant figures in net 
energy generation per unit mass and generally 2-3 significant figures in the mass fractions of 
various isotopes. Hence, the 8 isotope network is sufficiently accurate to describe the energy 
generation for hydrogen and helium burning. 



5. Macrophysics 
5.1. The mixing length theory of convection 

The mlt module implements the standard mixing length theory (MLT) of convection 



as presented by Cox & Giuli (1968, chapter 14). There are options both for computing the 



actual temperature gradient, Vt, when the total luminosity, L, is specified and for computing 
the convective luminosity, L conv , when Vt is specified. The mlt module calculates diffusion 
coefficients for those codes, such as MESA star, that treat convective mixing of elements as 
a diffusive process. The quantities listed in Tableland their partial derivatives with respect 



http : / / groups . nscl . msu . edu/ j ina/ reaclib/ db/ index . php 
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to several physical variables are returned by the mlt module. 



In addition to the standard MLT of Cox & Giuli, the mlt module includes the option 



to use the modified MLT of Henyey et al. (1965). Whereas the standard MLT assumes 



high optical depths and no radiative losses, the Henyey et al. (1965) variation allows the 



convective efficiency to vary with the opaqueness of the convective element, an important 



effect for convective zones near the outer layers of stars. If the Henyey et al. (1965) option 



is used, the parameter v (a mixing length velocity multiplier) and y (a parameter that sets 
the temperature gradient in a rising bubble) may be set by the user. They default to the 
recommended values of y — 1/3 and v = 8. 



Towards the center of a star, the commonly used definition of the pressure scale height, 
Xp = P/gp, diverges as g — > 0. Therefore, we provide the option of using the alternate 
definition of |Eggleton| ( |l97l[ ), X' P = (P/Gp 2 ) 1 / 2 , 
X' P ~ R. 



when X' P < Xp. 



At the center of the star, 
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Table 6. MESA mlt output quantities and units 



Output 


Definition 


Units 


V T 


actual temperature gradient a 


Dimensionless 


V r ad 


radiative temperature gradient 21 


Dimensionless 


-^conv 


convective luminosity 15 


ergs s _1 


L 


total luminosity 


ergs s _1 


Ap 


pressure scale height 


cm 


A 


mixing length (= cimlt^p) 


cm 


^conv 


convective velocity 


cm s _1 


D 


Eulerian diffusion coefficient 


2 —1 

cnr s 


a 


Lagrangian diffusion coefficient (= D(Airr 2 p) 2 ) 


g 2 s- 1 



a Only when L is specified. 

b Only when Vt is specified. 

c Only when Vt is specified and L conv > 0. 
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5.2. Convective overshoot mixing 



As described in §6.2[ MESA star treats convective mixing as a time-dependent, diffusive 
process with a diffusion coefficient, D, determined by the mlt module. In the absence of a 3-D 
hydrodynamical treatment of convection it is necessary to account for the hydrodynamical 
mixing instabilities at convective boundaries, termed overshoot mixing, via a parametric 
model. After the MLT calculations have been performed, MESA star sets the overshoot 
mixing diffusion coefficient 

Aw = Axmv.oexp ( —j^— J > ( 2 ) 
V / a p,o J 

where -D CO nv,o is the MLT derived diffusion coefficient at a user-defined location near the 
Schwarzschild boundary, X P0 is the pressure scale height at that location, z is the distance 



in the radiative layer away from that location, and / is an adjustable parameter (Herwig 



2000). In MESA star the adjustable parameter, /, may have different values at the upper 
and lower convective boundaries for non-burning, H-burning, He-burning, and metal-burning 
convection zones. 

Parameters are provided to allow the user to set a lower limit on -Dov below which 
overshoot mixing is neglected and to limit the region of the star over which overshoot mix- 
ing will be considered. So as to model the 13 C pocket needed for s-process nucleosynthesis, 
MESA star also allows an increase in the overshooting parameter at the bottom of the con- 



vective envelope during the third dredge-up compared to the inter-pulse value (Lugaro et al. 



2003). There is also an option to change the value of overshoot mixing at the bottom of the 
AGB thermal pulse-driven convection zone compared to the standard value chosen for the 
bottom of the He-burning convection zone. 



5.3. Atmosphere boundary conditions 



As described in §6.2[ the pressure, P s , and temperature, T s , at the top of the outermost 
cell in MESA star must be set by an atmospheric model. This is done by the atm module, 
which uses M, R, and L to provide P s and T s . It also gives partial derivatives of T s and P s 
with respect to the input variables. The atm module assumes the plane parallel limit, so that 
the relevant variables are g = GM/R 2 and T e g 4 = L/47to"sb-R 2 - With some options, the user 
must specify the optical depth r s to the base of the atmosphere, whereas in other cases, the 
atm module has an implicit value. Three methods are supplied by atm: direct integrations, 
interpolations in model atmosphere tables, and a "simple" recipe. 

The integrations of the hydrostatic balance equation, dP gabS /dr = g/n — (a/3)dT 4 /dr, 
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with dr = —npdr are performed using either the relation T 4 (r) 



3T e 4 ff (r+2/3)/4 (Eddington 



T 



1926), or the specific T — r relation of Krishna Swamy| (1966). These integrations start at 
10~ 5 and end at a user specified stopping point, r s , which defaults to t s = 2/3 (0.312) 
for Eddington (Krishna Swamy)] 1 "^] The routine integrates the gas pressure and then adds 
the radiation pressure at the stopping point to get P s . 

The MESA model atmosphere tables come in two forms. The MESA photospheric tables 
(which return T s = T c g and assume that r s w 1) cover log Z/Z Q = —4 to +0.5 assuming 
the Grevesse & Noels (1993) Solar abundance mixture. They span log(g) = —0.5 to 5.5 
at 0.5 dex intervals and T c g =2, 000-50, 000K at 250K intervals. They are constructed, in 



precedence order, with, first, the PHOENIX (Hauschildt et al. 1999a[b ) model atmospheres 
(which span log(g) = —0.5 to 5.5 and T cff = 2,000 to 10,000 K); and second, the Castelli 
& Kurucz (2003) model atmospheres (which span log(g) = to 5 and T e g = 3500 to 50, 000 
K). In regions where neither table is available, we generate the MESA table entry using the 
integrations described above with the Eddington T-r relation. The second MESA table is 
for Solar metallicity and gives P s and T s at r s = 100. It is primarily for the evolution of 
low mass stars, brown dwarfs, and giant planets. It is constructed from Castelli & Kurucz 



(2003), and for T eS < 3000K, the COND model atmospheres (Allard et al. 2001) which 



assume gravitational settling of those elements that form dust, depleting those elements 
from the photosphere. Figure [6] shows the regions where the different sources are used, and 
in those regions where there are no published results, we use the integration of the Eddington 
T-r relation. 



18 If the first attempt to integrate fails, the code makes two further attempts, each time increasing the 
initial r by a factor of 10. The integration is carried out with the Dormand-Price integrator from the num 
module. 
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Fig. 6. — The range of T e g and log(g) covered by the MESA atm tables for r s = 100 and Solar 
metallicity. The CK region uses the tables of |Castelli fc Kuruczj ( |20"03| , whereas the COND 
region uses Allard et al. (2001). At lower log(g) and cold regions, we use direct integrations 



of the Eddington T — t relation. The green lines show evolutionary tracks of stars, brown 
dwarfs and giant planets of the noted masses. 



Finally, there is a simple option where the user specifies r s and we use the constant 
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opacity, k s , solution of radiative diffusion, 



Ts9 



1 + 1.6 x 1(T 4 k s ( ^/ L ® 

\M/M G 



(3) 



where the factor in square brackets accounts for the nonzero radiation pressure (see, e.g., 
Cox fc Giuli||1968 Section 20.1). The temperature is simply given by the Eddington relation. 
The user can either specify k s or it will be calculated in an iterative manner using the initial 
value of P s from an initial guess at k s (usually given by MESA star as the value in the 



outermost cell; see £6.2). In addition, the atm module has the option to revert to Equation 



(J3| if a model wanders outside the range of the currently used model atmosphere tables or 
if the atmosphere integration fails for any reason. 



5.4. Diffusion and gravitational settling 



MESA diffusion calculates particle diffusion and gravitational settling by solving Burger's 



equations using the method and diffusion coefficients of Thoul et al. (1994). The transport 



of material is computed using the semi- implicit, finite difference scheme described by lb en &: 



MacDonald (1985). Radiative levitation is not presently included. The diffusion module 



treats the elements present in the stellar model as belonging to "classes" defined by the 
user in terms of ranges of atomic masses. For each class, the user specifies a representative 
isotope, and all members of that class are treated identically with their diffusion velocities 
determined by the representative isotope, and the diffusion equation solved with the mass 
fraction in that class. The caller can either specify the ionic charge for each class at each 
cell in the model or have the charge calculated by the ionization module, which estimates 
the typical ionic charge as a function of T, p, and free electrons per nucleon from |Paquette 



et al. (1986). 
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Fig. 7. — The absolute values of the diffusion velocities from diffusion (lines) and those 



published by Thoul et al. (1994). All results are plotted in units of Rq/tq, where tq = 
6 x 10 13 yr is the characteristic diffusion timescale for the Sun (Thoul et al. 1994). The 



dark solid and dashed lines are the diffusion results for H and O. The filled green circles 
show the results of Thoul et al. (1994) for H, O and Fe (Z = 21). The diffusion results 
for Helium are shown as the dashed red line. The diffusion results for Fe include one 
for Z = 21 (dotted blue line) and one for ionization states determined by ionization (the 
dot-dashed blue line). 
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The lines in Figure [7] plot four classes (H, He, 0, and Fe) with a solar model from 



MESA star and compares where possible to the results from Figure 9 of Thoul et al. (1994) 



shown by the filled green circles. The agreement is excellent for H, O and Fe (when we fix Fe 
to have the Z = 21 ionization state chosen by Thoul et al. (1994)). |Thoul et al.| ( |1994 ) did 
not exhibit the He velocity, so we have no comparison. For Fe, we also show the diffusion 
velocity when ionization finds a changing ionization state in the Z = 16,17,18 region 
(shown by the upper dot-dashed blue line), highlighting the need to better determine the Fe 



ionization state (Gorshkov & Baturin 2008). We also compared the diffusion output to 



the recent calculations of Gorshkov & Baturin (2008), finding agreement at better than 5% 
for the Fe case at Z = 26 and for O. 

The diffusion calculation can be restricted to areas where the temperature is above some 
minimum value, or where the mass fraction of a diffusing element is above some minimum 
value, aiding the convergence of solutions in a variety of environments. The physics imple- 
mentation is presently limited to regions where the Coulomb coupling parameter, T, is less 
than unity. At present, this inhibits an accurate calculation for segregation and settling of 
the remaining envelope H and He envelope on a cooling white dwarf. 



5.5. Testing MESA modules in an existing stellar evolution code 

The complex, nonlinear behavior of stellar structure and evolution models makes it 
difficult to disentangle the effects of model components (e.g., EOS, opacities, boundary 
conditions, etc.) when comparing results of separate codes. By design, the modularity of 
MESA allows individual physics modules to be incorporated into an existing stellar evolution 
code, tested, and then compared against the prior implementation of comparable physics in 
the same code. 

During the development of MESA, several MESA modules were integrated into the Dart- 
mouth Stellar Evolution Program (DSEP, Potter et al.|2007 ). This section reports the results 



of using four MESA modules, eos, kap, atm, and mlt, in DSEP to compute the evolution of 
a l.OM star with initial values of X = 0.70 and Z = 0.02. The star was evolved from the 
fully convective pre-main sequence to the onset of the core He flash. This was done six times: 
once, as the control case, using only DSEP routines and no MESA modules; next, using each 
of four MESA modules individually; and, finally, using the four MESA modules at the same 
time in DSEP. 

DSEP employs a p(P,T) EOS and so the MESA P gas — T tables were used during the 
eos test. Though DSEP and kap use the same sources for radiative opacities, they differ 
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in interpolation methods and the treatment of electron conduction opacities (see Bjork & 



Chaboyer 2006, for a thorough list of the physics in DSEP). When atm was tested, we used 



the Eddington grey atmosphere model integrated to r = 2/3. DSEP uses the Henyey et al. 



(1965) modification of the mixing length theory, which is available in mlt, and assumes that 
convective regions are instantaneously mixed to a uniform composition. 




Fig. 8. — Comparison of DSEP tracks using built-in physics modules and MESA modules 
for opacities, EOS, mixing length theory, and the atmospheric boundary condition. These 
tracks are for a l.OM star with initial X = 0.70 and Z = 0.02 evolved from the fully 
convective pre-main sequence to the onset of the He core flash. Only the H-R diagram shows 
the full evolutionary track. The T c panels omit the pre-main sequence in order to highlight 
the regions where the differences are most pronounced; the lifetime panel focuses on the end 
of the main sequence and red giant phase for the same reason. 

DSEP tracks employing either the atm or the mlt modules produce results that agree 
with the DSEP-only track to about 1 part in 10 4 . DSEP tracks employing the kap and eos 
modules exhibit some difference when compared to the DSEP-only track but, even in these 
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cases, the main sequence lifetime differs by less than 0.3% and T e g differs by less than 10K 
along the main sequence. As shown in Figure [8j the largest discrepancy between the DSEP- 
only track and the one that employs all four MESA modules appears in the T c — p c diagram 
when p c > 3 x 10 4 g cm -3 , corresponding to the growing helium core in the center of the red 
giant. Above logp c = 4, the track employing MESA modules is hotter than the DSEP-only 
track by ~ 0.02 in logT c at constant logp c . The center of the model has entered the region 
of electron degeneracy and electron conduction has become an important source of opacity. 
The majority of the difference is due to the EOS whereas the opacity difference amounts to 
about —0.005 in logT c , in the opposite direction to the EOS. The hotter conditions produced 
by the eos module is likely the cause for the slightly shorter RGB lifetime that can be seen 
in Figure |8j 



Stellar structure and evolution 



MESA star is a full- featured stellar structure and evolution library that utilizes the nu- 
merics and physics modules described in §'s [3]j5l It provides a clean-sheet implementation 



of a Henyey style code (Henyey et al. 1959) with automatic mesh refinement, analytic Ja- 



cobians, and coupled solution of the structure and composition equations. The design and 
implementation of MESA star was influenced by a number stellar evolution and hydrody- 



namic codes that were made available to us: EV (Eggleton 1971), EVOL (Herwig 2004), 



EZ (Paxton 2004), FLASH-the-tortoise (Lesaffre et al. 2006), GARSTEC (Weiss & Schlattl 



2008), NOVA (Starrfield et al. 2000), TITAN (Gehmeyr & Mihalas 1994), and TYCHO 



(Young & Arnett 2005). 



We now briefly describe the primary components of MESA star. MESA star first reads 



the input files and initializes the physics modules (see £6.1) to create a nuclear reaction 



network and access the EOS and opacity data. The specified starting model or pre-main 



sequence model is then loaded into memory (see £6.1), and the evolution loop is entered 



The procedure for one timestep has four basic elements. First, it prepares to take a new 



timestep by remeshing the model if necessary (£6.5 and 6.4). Second, it adjusts the model 
to reflect mass loss by winds or mass gain from accretion (£6.6) , adjusts abundances for 



element diffusion (£ |5.4 ), determines the convective diffusion coefficients (£ |5.1 and 5.2), and 



solves for the new structure and composition (£6.2 and 6.3) using the Newton- Raphson solver 
J3~l). Third, the next timestep is estimated (£6.4). Fourth, output files are generated (£6.1). 
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6.1. Starting models and basic input/output 

MESA star receives basic input from two Fortran namelist files. One file specifies the 
type of evolutionary calculation to be performed, the type of input model to use, the source of 
EOS and opacity data, the chemical composition and nuclear network, and other properties 
of the input model. The second file specifies the controls and options to be applied during 
the evolution. 

There are two ways to start a new evolutionary sequence with MESA star. The first 
is to use a saved model from a previous run. A variety of saved models are distributed 
with MESA as a convenience. These saved models fall into three general categories: (1) Zero 
Age Main Sequence (ZAMS) models for Z = 0.02 with 32 masses between 0.08 and 1OOM 
(MESA star will automatically interpolate any mass within this range); (2) very low mass, 
pre-main sequence models for Z = 0.02 and masses from 0.001 to O.O25M ; and (3) white 
dwarf models for Z = 0.02 with He cores of 0.15 - O.45M , C/O cores of 0.496 - l.O25M , 
and O/Ne cores of 1.259 — 1.376M . The user can also create saved models for essentially 
any purpose through available controls. 

The second way to start a new evolution is to create a pre-main sequence (PMS) model 
by specifying the mass, M, a uniform composition, a luminosity, and a central temperature, 
T c low enough that nuclear burning is inconsequential (T c = 9 x 10 5 K by default). For a fixed 
T c and composition, the total mass depends only on the central density, p c . An initial guess 
for p c is made by using the n — 1.5 polytrope, which is appropriate for a fully convective 
star, but we do not assume the star is fully convective during the subsequent search for 
a converged PMS model. Instead, MESA star uses the mlt, eos, and Newton solver from 
num to search for a p c that gives a model of the desired mass. The PMS routine presently 
creates starting models for 0.02 < M/M & < 50. Beyond these limits we find challenges 
converging the generated PMS model within the MESA star evolutionary loop. For such 
cases it is currently better to generate a starting model within the acceptable mass range, 



save it, relax it to a new mass with a specified mass gain or loss (see £6.6), and save that 
model. 

MESA star has the ability to create a binary file of its complete current state, called a 
photo, at user-specified timestep intervals. Restarting from a photo ensures no differences 
in the ensuing evolution. When restarting from a photo, many controls and options can 
be changed. A photo is different than a saved model in that a saved model is a text file 
containing a minimal description of the structure and composition but does not have enough 
information to allow a perfect restart. However, saved models are not tied to a particular 
version of the code and therefore are suitable for long term use or sharing with other users. 
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There are two additional types of output files, logs and profiles. A log records evolu- 
tionary properties over time such as stellar age, current mass, and a wide array of other 
quantities. A profile records model properties at a specified timestep at each zone from sur- 
face to center. MESA star can also output models in the FGONG formal 19 for use with stellar 



pulsation codes and se output for nucleosynthesis post-processing with NuGrid codes ^ Fi- 
nally, a few simple lines of user-supplied code allows for saving variables or combinations of 
variables that are not in the list of predefined options. 



6.2. Structure and composition equations 

MESA star builds 1-D, spherically-symmetric models by dividing the structure into cells, 
anywhere from hundreds to thousands depending on the complexity of nuclear burning, gra- 
dients of state variables, composition, and various tolerances. Cells are numbered starting 
with one at the surface and increasing inward. MESA star does not require the structure 
equations to be solved separately from the composition equations (operator splitting). In- 
stead, it simultaneously solves the full set of coupled equations for all cells from the surface 
to the center. The solution of the equations is done by the Newton solver from num using 
either banded or sparse matrix routines from mtx. The partial derivatives for use by the 
solver are calculated analytically using the partials returned by modules such as eos, kap, 
and net. 



|http : // owww . phys . au . dk/~ j cd/solar_mo dels/ 
2r http : //forum . astro . keele .ac.uk: 8080/nugrid 
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Fig. 9. — Schematic of some cell and face variables for MESA star. 



Each cell has some variables that are mass-averaged and others that are defined at 
the outer face, as shown in Figure [9j This way of defining the variables is a consequence 
of the finite volume, flux conservation formulation of the equations and improves stability 



and efficiency (Sugimoto et al. 1981). The inner boundary of the innermost cell is usually 



the center of the star and, therefore, has radius, luminosity, and velocity equal to zero. 
Nonzero center values can be used for applications that remove the underlying star (e.g., the 
envelope of a neutron star), in which case the user must define the values of M c and L c at 
the inner radius R c . The cell mass-averaged variables are density p k , temperature T k , and 
mass fraction vector X itk . The boundary variables are mass interior to the face m k , radius 
r k , luminosity L k , and velocity v k . In addition to these basic variables, composite variables 
are calculated for every cell and face, such as e nuc , k, a k , and F k (see Table [l] for variable 
definitions). All variables are evaluated at time t + 5 t unless otherwise specified. 
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The density evolution of cell k is determined by a finite volume form of the mass con- 
servation equation 



Pk 



dm i 



(4/3)7r(r2 - r^j) 



(4) 



For the innermost cell, r^+i is replaced by the inner boundary condition which is typically 
zero but can be nonzero for some applications. We reformulate many of our equations to 
improve numerical stability of the linear algebra and minimize round-off errors. We thus 
rewrite equation Q as 

1 , f 3 3 drrik 



log r k = - log 



fe+i 



4vr p k 



(5) 



The velocity of face k is zero unless the hydrodynamics option is activated, in which 



case 



d()ogr k ) 
dt 



(6) 



is the Lagrangian time derivative of the radius at face k. For enhanced numerical stability, 
we rescale this equation by dividing by the local sound speed. 

The pressure Pk is set by momentum conservation at interior cell boundaries, 



Pk-l — Pk 



drrik 

dlTLk 



(-) 

\dmj hydrostatic 



\ / hydrodynamic 



Aixrl 



(7) 



where drrik — 0.5(drrik-i + drrik), and a& is the Lagrangian acceleration at face k, evaluated 
by the change in v k over the timestep St. The acceleration is set to zero if the hydrodynamic 
option is not used. Similarly, the temperature of interior cells Tk is set by energy transport 
across interior cell boundaries, 



-k-x 



drrii 



\dmj hydrostatic Pk 



where Vt,/c = d log T/d log P at face k from the MESA module mlt (see { 5.1 ), Tk = (Tk-idrrik + 
Tkdrrik-i) / (drrik + dm^-i) is the temperature interpolated by mass at face k, and Pk = 
(Pk-idrrik + Pkdrrik-i) / (drrik + d> m k-i) is the pressure interpolated by mass at face k. For 



enhanced numerical stability, we rescale equation (|7j) by dividing by Pf, and equation (|8j) by 
dividing by T k - 

The pressure and temperature boundary conditions are constructed by using P s and T s 



from the MESA module atm (see £5.3). The difference in pressure and temperature from the 
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surface to the center of the first cell is found from hydrostatic equilibrium and Vt by 

Gmidmi/2 



dP x 



i 



dT s = dP s V T ,i^ ■ (9) 
n 

The boundary conditions are then 

log 71 = \og(T s + dT s ) 

log Pi = \og(P s + dP s ) . (10) 

These implicit equations for Pi and T\ are solved together with the regular structure and 
composition equations. 

Our finite volume form of energy conservation for cell k is 

Lk ~ Lk+i = dl7lk{€ nn c ~ Ci/,thermal + Cgrav) j (11) 

where e nuc (from module net or j ina) is the total nuclear reaction specific energy generation 
rate minus the nuclear reaction neutrino loss rate, and e^thermai (from module neu) is the 
specific thermal neutrino loss rate. The e grav term is the specific rate of change of gravitational 
energy due to contraction or expansion, 



, rflogT dlogp 



(12) 



where dlogT/dt and dlogp/dt are Lagrangian time derivatives at cell center by mass, and 
the other symbols are defined in Tables [3] and [6j For the innermost cell, L^+i is replaced 
by the inner boundary condition which is typically zero but can be nonzero, L c , in specific 
applications. For additional numerical stability, we rescale equation (11) by dividing by a 
scale factor that is typically the surface luminosity of the previous model. 

The equation for mass fraction of species i in cell k is 
Xi 7 k(t + St) — Xi^it) = dX hurn + dX m i x 

where dXi^/dt is the rate of change from nuclear reactions reported by net or jina, is 
the mass of species i flowing across face k 

Fi,k = {Xi^k — Xi t k-i) ==- , (14) 

arrik 
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where is the Lagrangian diffusion coefficient from the combined effects of convection 
(§ |5.1 ) and overshoot mixing (£5.2). For numerical stability, is calculated at the beginning 
of the timestep and held constant during the implicit solver iterations. This assumption 
accommodates the non-local overshooting algorithm and significantly improves the numer- 
ical convergence. It leads to a small inconsistency between the mixing boundary and the 
convection boundary as calculated at the end of the timestep. 

Equations ([5]), d7|), (|8j, (11), (14), and, optionally equation (|6]), are by default solved 
fully coupled and simultaneously with a 1 st order backwards differencing time integration. 



6.3. Convergence to a solution 



The generalized Newton-Raphson scheme is represented by 



= F{y) = F{y t + 6%) = F(yJ + 

where is a trial solution, F(y*i) is the residual, Sy, 
Jacobian matrix. 



dF 
dy 



Sy t + 0{8yl 



(15) 



is the correction, and [dF/dy\i is the 



MESA star uses the previous model, modified by remeshing, mass change, and element 
diffusion, as the initial trial solution for the Newton-Raphson solver. This is generally suc- 



cessful because we use analytic Jacobians and have sophisticated timestep controls (see £6.4). 
The use of analytic Jacobians in MESA star requires that each of the MESA modules provides 
not just the required output quantities but also quality, preferentially analytic, partial deriva- 
tives with respect to the input quantities. At each timestep, MESA star converges on a final 
solution by iteratively improving upon the trial solution. We calculate the residuals, con- 
struct a Jacobian matrix, and solve the resulting system of linear equations with the solvers 
in mtx to find the corrections to the variables. 

The trial solution is accepted when the corrections and residuals meet a specifiable set 
of comprehensive convergence criteria. In most cases, the solver is able to satisfy these limits 
in 2 or 3 iterations. However, under difficult circumstances like the He core flash or advanced 
nuclear burning in massive stars, MESA star can automatically adjust the convergence cri- 
teria. The corrections to the variables will, generally, not produce zero residuals because the 
system of equations is nonlinear. In some cases, the corrections might make the residuals 
larger rather than smaller. In such cases, the length of the correction vector is reduced by a 
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line search scheme, 2 *| until they improve the residuals. In principle, the residuals can be made 
arbitrarily small, but this may take a prohibitively large number of iterations. In practice, 
the use of the line search scheme helps the convergence rate in many cases, but cannot ensure 
convergence in all cases. 

If convergence cannot be achieved with the current timestep, then MESA star will first 
try again with a reduced timestep (a "retry" ) anticipating that a smaller timestep will reduce 
the non-linearity. If the retry fails, MESA star will return to the previous model and with 
a smaller timestep than it used to get to the current model (a "backup"). If the backup 
fails, MESA star will continue to reduce the timestep until either the model converges or 
the timestep reaches some pre-defined minimum, in which case the evolutionary sequence is 
terminated. 



6.4. Timestep selection 



Timestep selection is a crucial part of stellar evolution. The timesteps should be small 
enough to allow convergence in relatively few iterations, but large enough to allow efficient 
evolutions. Changes to the timestep should also provide for rapid responses to varying struc- 
ture or composition conditions, but need to be carefully controlled to avoid over-corrections 
that can reduce the convergence rate. 

MESA star does timestep selection as a two stage process. The first stage proposes 



a new timestep using a scheme based on digital control theory (Soderlind & Wang 2006). 



The second stage implements a wide range of tests that can reduce the proposed timestep 
if certain selected properties of the model are changing faster than specified. For the first 
stage, we use a low-pass filter. The control variable v c is the unweighted average over all cells 
of the relative changes in logp, logT, and \ogR. The target value v t is 10~ 4 by default. For 
improved stability and response, the low-pass filter method uses the previous two results. 
Let 8ti, and 5U+i be the previous, current, and next timestep, respectively, while w c ,i-i 

and v C) i are the previous and current values of v c . The timestep for model % + 1 is then 
determined by 



SU 



i+1 



SUf 



1/4 



(16) 



fidk/dti-x) 

where f(x) = 1 + 2tan _1 [0.5(x — 1)]. The control scheme implemented by equation (16) 
allows rapid changes in timestep without undesirable fluctuations. 



1 This is a globally convergent method and is similar to what is described in §9.7 of 



Press et al. 



(1992). 
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The timestep proposed by this low-pass filtering scheme can be reduced according to a 
variety of special tests that have hard and soft limits. If a change exceeds its specified hard 
limit, the current solution is rejected, and the code is forced to do a retry or a backup. If a 
change exceeds its specified soft limit, the next timestep is reduced proportionally. Examples 
of special tests include limits on the maximum absolute or relative changes in mesh structure, 
composition variables, nuclear burning rate, T eS , L, M, T c , p c , and integrated luminosity 
from various types of nuclear burning. 



6.5. Mesh adjustment 

MESA star checks the structure and composition profiles of the model at the beginning of 
each timestep and, if necessary, adjusts the mesh. Cells may be split into two or more pieces, 
or they may be made larger by merging two or more adjacent cells. The overall remeshing 
algorithm is designed such that most cells are not changed during a typical remesh. This 
minimizes numerical diffusion and tends to help convergence. Remeshing is divided into a 
planning stage and an adjustment stage. 

The planning stage determines which cells to split or merge based on allowed changes 
between adjacent cells. Mesh revisions minimize the number of splits and maximize the 
number of merges while ensuring that the magnitudes, A, of differences between any two 
adjacent cells are below specific thresholds: AlogP < Op, AlogT < T , and A log[X( 4 He) + 
X( 4 Heo)] < #He where X( 4 He) is the helium mass fraction and X( 4 Heo) sets an effective lower 
lower limit on the sensitivity to the helium abundance. The default thresholds are Op = 1/30, 
0t = 1/80, #He = 1/20, and X( 4 Heo) = 0.01. Options are available for specifying allowed 
changes between cells for other mass fractions, AV n d and Alog(T/(T + T )) for arbitrary 

Local reductions in the magnitude of allowed changes will place higher resolution in 
desired regions of the star. For example, the default is to increase resolution in regions of 
nuclear burning having A log e nuc large compared to A log P. This increase takes effect at a 
minimum loge nuc = —2 and increases to a maximum factor of 4 in resolution for loge nuc > 4. 
The size and range of enhancement can also be set for various specific types of burning. 
Similarly, it is possible to increase resolution near the boundaries of convection zones over a 
distance measured in units of the pressure scale height. Different enhancements and distances 
can be specified for above and below the upper and lower boundaries of zones with or without 
burning. There are also options to increase spatial resolution in regions having A log Xi large 
compared to A log P, or near locations where there are spatial gradients in the most abundant 
species. Finally, further splitting is done as necessary to limit the relative sizes of adjacent 
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cells. 

The adjustment stage executes the remesh plan. Cells to be split are constructed by first 



performing a monotonicity preserving cubic interpolation (Steffen 1990) in mass to obtain 
the luminosities and enclosed volumes at the new cell boundaries. The new densities are 
then calculated from the new cell masses and volumes, as shown in equation Q. Next, new 
composition mass fraction vectors are calculated. For cells being merged, this is straightfor- 
ward. For cells being split, neighboring cells are used to form a linear approximation of mass 
fraction for each species as a function of mass coordinate within the cell. The slopes are 
adjusted so that the mass fractions sum to one everywhere, and the functions are integrated 
over the new cell mass to determine the abundances. 

Finally, the method for calculating the new temperature varies according to electron 
degeneracy. As the electrons become degenerate (i.e. 77 > 0), split cells simply inherit their 
temperature while merged cells take on the mass-average of their constituent temperatures. 
If the electrons are not degenerate (i.e. rj < 0), then a reconstruction parabola is created 



for the specific internal energy profile of the parent and its neighbor cells (Stiriba 2003). 
The parabola is integrated over the new cell to find its total internal energy. The new cell 
temperature is determined by repeatedly calling the eos module using the new composition 
and density with trial temperatures until the desired internal energy is found. 



6.6. Mass loss and accretion 



Mass adjustment for mass loss or accretion is done at each timestep before solving the 
equations for stellar structure and composition. MESA star offers a variety of ways to set 
the rate of mass change M. A constant mass accretion or mass loss rate may be specified in 



the input files (see {6.1). Implementations of Reimers (1975) for red giants, Blocker (1995) 



for AGB stars, de Jager et al. (1988) for a range of stars in the H-R diagram, mass loss for 



massive stars by (Glebbeek et al. 2009; Vink, de Koter & Lamers 2001 Nugis & Lamers 



2009[ |Nieuwenhuijzen fc de Jager||1990~l ), supersonic mass loss inspired by iPrialnik & Kovetz 
(1995), and super-Eddington mass loss ( Paczynski fc Proszynski|[l986 ) are available options. 
An arbitrary mass accretion or mass loss scheme may be implemented by writing a new 



module. An example of such a routine provided with MESA star is Mattsson et al. (2010) 



mass loss for carbon stars. Finally, one may write a wrapper program that calculates M for 
each timestep and then calls the MESA star module. 

Since MESA star allows for simulations with a fixed (and unmodeled) inner mass, M c , 
the total mass is M = M c + M m , where M m is the modeled mass. For cell k, MESA star 
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stores the relative cell mass dq k = dm k /M m and the relative mass interior to a cell face 



Qk = m k/M m = 1 — Yll=i 1 dqi (see Figure 9). Rather than evaluate dq k as q k — qk+i, it is 



essential to define q in terms of dq to maintain accuracy (Lesaffre et al. 2006). For example, 
in the outer envelope of a star where the q k approach 1, the dq k can be 10 -12 or smaller. 
Subtraction of two adjacent q k to find a dq k leads to a intolerable loss of precision. 

After a change in mass, 5M, has been determined, the mass structure of the stellar 
model is modified. This procedure changes the mass location of some cells and revises the 
composition of those cells to match their new location. It does not add or remove cells, 
nor does it change the initial trial solution for the structure variables such as p, T, r, or 
L. The mass structure is divided into an inner (usually the central regions of the star), an 
intermediate, and an outer region (usually the stellar envelope). The boundaries of the inner 
and outer regions are initially set according to temperature, with defaults of logT = 6 for the 
inner boundary and log T = 5 for the outer boundary. This range is automatically expanded, 
for enhanced numerical stability, if the mass in the intermediate region is not significantly 
larger than 5M. The range is first enlarged by moving the outer boundary to the surface. If 
the enclosed mass in the intermediate region is still too small, then the inner boundary can 
be moved inward subject to certain limits. One limit is that the inner boundary does not 
cross a region of the model where the composition changes rapidly. Another limit is that the 
fractional mass of the intermediate region cannot change by more than a factor of two from 
its previous value nor exceed 10% of the total mass. 

Once the regions have been defined, the dq k are updated. In the inner region the dq k 
are rescaled by MjiM + SM). Thus, dm k , m k , and X k have the sames value before and after 
a change in mass to eliminate the possibility of unwanted numerical mixing in the center. 
In the outer region, cells retain the same value of dq k to improve convergence in the high 



entropy regions of the star (Sugimoto et al. 1981). The dq k in the intermediate region are 



scaled so that ^2 dqk = 1- The composition of cells in the intermediate and outer regions are 
then updated. In the case of mass accretion, the composition of the outermost cells whose 
enclosed mass totals SM is set to match the specified accretion abundances. Cells that were 
part of the old structure have their compositions set to match the previous composition. 



6.7. Resolution sensitivity 

We examined the resolution convergence properties of a 1M model by varying the 
parameters for mesh refinement and timestepping. The mesh refinement parameter multi- 
plies the limits for variable changes across mesh cells and is closely correlated with the cell 
size. The timestepping parameter controls the tolerance of the cell average of the relative 
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changes between time steps in logp, logT, and \ogR (see £6.4) and is closely correlated 
with the timestep. We varied the mesh refinement and timestepping controls in tandem 
through a parameter C, which is a multiplicative factor on their default values of 1 and 10~ 4 , 
respectively. Therefore, C is anti-correlated with the time and space resolution. 

Table [7] and Figure 10 detail the convergence properties with C of a solar metallicity, 



1.0 M Q model with an 7^=0.5 Reimers mass loss model (Reimers 1975, see £6.6). These 
calculations begin at the ZAMS and are terminated at 11.0 Gyr, when the model stars 
are turning off the main sequence. As a measure of convergence, we use the difference, £, 
between a quantity at a given resolution and the quantity at the highest resolution considered 
(C=l/16). In order to determine how convergence depends on resolution (|£| oc C°), we 
determine the order of convergence, a, for increasingly resolved pairs in Table [7} 

a = log(p^) /\og(-^) . (17) 



^coarse / ' \ C 

The convergence orders show that all values converge linearly at large values of C and 
display super-linear convergence (a ~ 1.6) at smaller values of C. These convergence orders 
are plausible given that we use a first order time integration scheme and a finite volume 
differencing scheme that is second order accurate on uniform grids. 



Table [8] and Figure 11 detail the same stellar models as a function of C except the 



calculations are stopped at L = 1OOL , when the stars are on the RGB. Table [8] suggests 
the age of the star converges linearly at larger values of C and super-linearly (a ~ 1.5) at 
smaller values of C. However, T e g, T c , p c , M, and Mn e all display oscillatory behavior about 
the C=l/16 solution, suggesting factors other than spacetime resolution are dominating the 
error at this stage of the evolution. Such factors could be limits in the precision attained 
by interpolation in the various tables (e.g., 4 significant figures for opacities, ~ 6 significant 
figures for the OPAL and SCVH EOS), or small changes in boundary conditions. 



The lower panel of Figure [12] shows the sound speed profile for the 100L Q model with 
C=l/16, our highest resolution case. The helium core and convective zone boundaries are 
labeled. The mass interior to the convective zone boundary is 0.95 M Q . The upper panel of 
Figure [12] shows the convergence properties of the sound speed profile with resolution in the 
hydrogen layer. Both the C=l and C=l/2 profiles have sound speeds that are smaller than 
the C=l/16 profile, while the C=l/4 and C=l/8 profiles have larger sound speeds. Note 
that the difference between the various profiles becomes less as the parameter C is made 
smaller, indicating that the convergence rate is becoming smaller. This suggests factors 
other than spacetime resolution are dominating the convergence rates at this stage of the 
evolution. Again, this could be due to the precision attained by table interpolations, small 
changes in boundary conditions, or differencing errors. 



-48 - 



Table 7. 1M Model Convergence Properties at 11.00 Gyr 



Control parameter C 


2 


1 


1/2 


1/4 


1/8 


1/16 


Number of cells 


457 


732 


1385 


2740 


5426 


10777 


Number of timesteps 


93 


135 


225 


418 


813 


1608 


L (£©) 


2.06094 


2.04251 


2.03241 


2.02678 


2.023737 


2.02217 


e (xio- 3 ) 


19.17 


10.06 


5.06 


2.28 


0.775 


0.0 


a 




0.93 


0.99 


1.15 


1.56 




T cff (K) 


5543.195 


5573.935 


5587.434 


5593.334 


5596.064 


5597.361 


€ (xio- 3 ) 


-9.68 


-4.19 


-1.77 


-0.72 


-0.232 


0.0 


a 




1.21 


1.24 


1.30 


1.63 




LogT c 


7.301645 


7.298305 


7.296797 


7.296052 


7.295676 


7.295486 


£ (xlO- 3 ) 


0.8442 


0.3864 


0.1797 


0.0776 


0.0260 


0.0 


a 




1.13 


1.10 


1.21 


1.57 




Log /9 C 


3.393658 


3.348505 


3.32615 


3.314372 


3.308441 


3.305552 


£ (xio- 3 ) 


26.65 


12.99 


6.23 


2.69 


0.874 


0.0 






1.04 


1.06 


1.22 


1.61 
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Fig. 10. — Convergence, £, for L, T e g, T c , and p c for a 1M model at 11.0 Gyr as a function 
of the control parameter C. These differences are all with respect to the C=l/16 model. 
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Fig. 11. — Convergence properties in stellar age, T e s, T c , p c , M, and Mn e for a Mi = 1M 
model at 1OOL as a function of the control parameter C. These differences are all relative 
to the C=l/16 model. Factors other than the spacetime resolution are dominating the 
differences for quantities other than the stellar age. 
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Table 8. 1M Q Model Convergence Properties at 100L Q 



Control parameter C 


Z 


1 
1 


1 /i 
1/2 


1 l A 

1/4 


1 /c 
1/8 


1 /1 c\ 

1/ lb 


Number of cells 


763 


1616 


3262 


6550 


13146 


26248 


Number of timesteps 


689 


1181 


2291 


4547 


8992 


17812 


Age (Gvrl 


12.302 


12.367 


12.400 


12.419 


12.428 


12.433 


€ (xicr 3 ) 


-10.54 


-5.31 


-2.65 


-1.13 


-0.402 


0.0 






0.99 


1.00 


1.24 


1.49 




Tpff (K) 

en v-*-*-/ 1 


4173.850 


4183.587 


4185.450 


4184.537 


4184.823 


4185.260 


e (xio- 3 ) 


-2.73 


-0.40 


0.05 


-0.17 


-0.10 


0.0 


a 




2.77 


3.14 


-1.93 


0.73 




Log T c 


7.61573 


7.61327 


7.61421 


7.613652 


7.613793 


7.613556 


i (xlO- 3 ) 


0.29 


-0.04 


0.09 


0.01 


0.03 


0.0 


a 




2.93 


-1.19 


2.77 


-1.30 




Log p c 


5.42340 


5.42393 


5.42402 


5.424929 


5.424868 


5.424713 


Z (XlO- 3 ) 


-0.24 


-0.14 


-0.13 


0.04 


0.03 


0.0 


a 




0.75 


0.18 


1.68 


0.48 




M (M ) 


0.984444 


0.984482 


0.984416 


0.984356 


0.984345 


0.984351 


e (xio- 3 ) 


0.09 


0.13 


0.07 


0.01 


-0.01 


0.0 


a 




-0.49 


1.01 


3.70 


-0.26 




M Hc (M ) 


0.28080 


0.28062 


0.28085 


0.281016 


0.281039 


0.281001 


€ (xio- 3 ) 


0.09 


0.13 


0.07 


0.01 


-0.01 


0.0 






-0.49 


1.01 


3.70 


-0.26 
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Fig. 12. — Sound speed profile (lower panel) and convergence properties (upper panel) for 
the 1M model at 1OOL with respect to the reference case, C=l/16. The helium core 
boundary and convective zone boundary are labeled. 



6.8. Multithreading 



MESA modules are designed to be thread-safe (see thereby enabling parallel execu- 
tion. Table [9] lists the execution times in seconds of several specific tasks from 4 identical 
evolutionary calculations, each with a different number of threads (one thread per core). The 
essential difference between modules that scale as the inverse of the number of threads (eos 
and net) and those that don't (e.g., kap and neu) is the ratio of the overhead associated 
with parallel execution to the actual time required for each module to perform its calculation. 
The time required to complete serial tasks sets a lower limit on the total execution time of 
MESA star. 



6.9. Visualization with PGstar 



By default, MESA star provides alpha- numeric output at regular intervals. In addition, 
it provides the option for concurrent graphical output. PGstar uses the PGPLOTQ library 
to create on-screen plots or images in PNG format that can be post-processed into anima- 
tions of an evolutionary sequence. A wide variety of visualization options are provided and 
these are all configurable through the PGstar inlist. For example, the PGstar window can 
simultaneously hold: an H-R diagram, a T c — p c diagram, and interior profiles of physical 
variables, such as nuclear energy generation, and composition. Animation is very useful for 
visualizing complex, time dependent processes. For example, view the short selection from 
the MESA website that shows the He core flash in a 1M star, 23 



7. MESA star results: comparisons and capabilities 

As with any modeling approach, MESA star must be verified ( "Is it solving the equations 
and validated ( "Does it solve the right equations?" ) to demonstrate 

. V&V is a maturing discipline 
, with the goal of assessing the error and uncertainty in 
a numerical simulation, which also includes addressing sources of error in theory, experiment, 

. The results of V&V testing are historical 
statements of reproducible evidence that a simulation demonstrates a quantified level of 



correctly?" Roache|1998 



its accuracy and predictive credibility (e.g., Oberkampf||2004 



e.g., Roache|1998 Calder et al.pOOl 



observation, and computation (Calder et al.||2004 



22 We thank Philip Pinto for initiating this MESA capability, see 



http: //www. astro . caltech. edu/~t jp/ 



pgplot/ 

23 



http: //mesa. sour ceforge .net /pdfs/lMHef lash. mov 
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Table 9. Execution times (s) with multiple threads 



Number of Threads 





1 


2 


4 


8 


total a 


12.2146 


8.0099 


5.8963 


5.0634 


ratio 




1.53 


1.36 


1.16 


Threaded tasks 


net 


6.2602 


3.1721 


1.6047 


0.8182 


eos 


1.7185 


0.8897 


0.4539 


0.2399 


mlt 


0.2479 


0.1384 


0.0704 


0.0357 


kap 


0.2285 


0.1386 


0.1044 


0.0875 


neu 


0.0240 


0.0209 


0.0139 


0.0098 


subtotal 


8.4791 


4.3597 


2.2473 


1.1911 


ratio 




1.95 


1.94 


1.89 


fraction of total 


0.69 


0.54 


0.38 


0.24 


Serial tasks 


file output 


1.1848 


1.0301 


1.0036 


1.1679 


matrix linear algebra 0.6569 


0.7654 


0.7885 


0.7926 


miscellaneous 


1.8938 


1.8547 


1.8569 


1.9118 


subtotal 


3.7355 


3.6502 


3.6490 


3.8723 


fraction of total 


0.31 


0.46 


0.62 


0.76 



a These numbers do not include initialization, e.g., loading 
of data tables. 
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accuracy in the solution of a specific problem. 



V&V is an ongoing activity for MESA via the MESA test suite (see Appendix |BJ), where 
code modules are tested individually, and, where possible, the integrated code MESA star 
is verified and validated. Verification for MESA includes a systematic study of the effect of 



mesh and time-step refinement on simulation accuracy (£6.7), specific module comparisons 
(§ |5.5[ ), and stellar evolution code comparisons presented in this section. 

This section shows MESA star evolution calculations of single stellar and substellar 
objects with 1O" 3 M < M < 1OOOM (in §§[7711 [7T2I O) as well as verification results 



£7.1.1, 7.1.2, 7.2.1, 7.2.2, 7.3.1, and 7.3.2). In £7.1.3 we compare the MESA star Solar 



model with helioseismic data. As examples of the many other experiments that are possible 
with MESA, we model prolonged accretion of He onto a neutron star and a mass-transfer 



scenario relevant to cataclysmic variables in £ 7.4 



7.1. Low mass stellar structure and evolution 

MESA star has sufficiently broad input physics to compute the evolution of low mass 
stars and substellar objects down to Jupiter's mass (~ 1O~ 3 M ), as well as complete evolu- 
tionary sequences of low mass stars (M < 2M ) from the PMS to the white dwarf cooling 



curve without any intervention. Figure 13 shows evolutionary tracks in the H-R diagram 



for 1 and 1.25M models with Z = O.Olrn The 1.25M m model exhibits a late He-shell flash 







during the pre- white dwarf phase. 



Figure 14 provides further examples, spanning O.9-2M at Z — 0.02; for clarity, the 
pre-main sequence portion of the tracks were removed and the runs were terminated after 
the models left the thermally-pulsating asymptotic giant branch (TP-AGB). The bottom 
panel shows the evolution in the T c — p c plane, exhibiting the convergence of the 0.9, 1.2 and 
1.5 M Q models to nearly identical, degenerate, Helium cores when on the RGB. The 2M 
model ignites He at a lower level of degeneracy. 



Each calculation takes a few hours on a laptop computer. 
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Fig. 13.— The MESA star evolution of a 1M and a 1.25M , Z = 0.01 star from the 
pre-main sequence to cooling white dwarfs. 
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Fig. 14.— Evolution from MESA star of 0.9, 1.2, 1.5, and 2 M stars with Z = 0.02 up to 
the end of the TP-AGB. The top panel shows their evolution in the H-R diagram, where the 
solid red point is the ZAMS. The bottom panel shows the evolution in the T c — p c plane, 
exhibiting the He core flash and later evolution of the C/O core during the thermal pulses. 
The dashed blue (heavy grey) line shows a constant electron degeneracy of e^/fc^T = 20(4). 
The dashed red line is for a constant pressure of logP = 20.25; relevant to the He core flash. 
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Though the required MESA star timesteps get short (~ hours) during the off-center He 
core flash in the 0.9, 1.2 and 1.5 M models, the stellar model does not become dynamic; 



the entropy change timescales are always longer than the local dynamical time (Thomas 



1967 Serenelli & Weiss 2005; Shen & Bildsten 2009 Mocak et al. 2009). The reduction of 



hydrostatic pressure in the core at the onset of flashes leads to the adiabatic expansion of 



the core, visible as the drop in T c at constant degeneracy. Successive He flashes (Thomas 
1967 Serenelli & Weiss 2005) work their way into the core over a 2 x 10 6 year timescale, 



eventually heating it (at nearly constant pressure; see the dashed red line in the bottom 



panel of Figure 14 ) to ignition and arrival onto the horizontal branch. The further evolution 



during He burning and the thermal pulses is seen in the bottom panel, where the small 
changes in the C/O core during thermal pulses on the AGB are resolved in MESA star. 



We start the more detailed calculations and comparisons to previous work in §7.1.1 by 
displaying the MESA star PMS evolution of low mass stars, brown dwarfs and giant planets, 



and comparing to prior results of (Baraffe et al. 2003). 



7.1.1. Low mass pre-main sequence stars, contracting brown dwarfs and giant planets 



The PMS evolution of low-mass stars (Burrows, Hubbard & Lunine 1989 D'Antona & 



Mazzitelli 


1994 


Baraffe et al. 


1998 


) and gravitationally contracting brown dwarfs and giant 


planets ( 


Burrows et al. 


1997 


Chabrier et al. 


2000 


Chabrier & Baraffe 


2000; 


Burrows et al. 



2001) has been studied extensively. The lasting importance of these problems motivates us 



to ensure that MESA star can successfully perform these evolutions. 



Figure 15 shows the evolution of PMS stars with masses of O.O8M < M < 1M . Each 
solid line starts on the PMS Hayashi track for a fixed mass, M, and ends at an age of 1 Gyr. 
All stars with M > O.2M have reached the ZAMS (L = L nuc ) by this time. The red circles 
show the location where the 7 Li is depleted by a factor of 100, but only for M < O.5M . Stars 



with M > O.5M deplete their 7 Li after the core becomes radiative (D'Antona & Mazzitelli 



1994 Chabrier, Baraffe & Plez 1996 Bildsten et al. 1997), adding an uncertain dependence 



on convective overshoot that we do not investigate here. 
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Fig. 15. — Location in the Hertzsprung-Russell (H-R) diagram for O.O85M < M < 1M & 
stars as they arrive at the main sequence for Y = 0.275 and Z = 0.019. The mass of the star 
is noted by the values at the bottom of the line. The dashed blue lines are isochrones for 
ages of 3, 10, 30, 100 and 300 Myr, as noted to the right. The purple squares (red circles) 
show where D ( 7 Li) is depleted by a factor of 100. The green triangles show the ZAMS. 

Lower-mass stars (M < O.3M ) remain fully convective throughout their PMS and 
arrival on the main sequence. In Figure [16] we show the evolution in the T c — p c plane 

1 /3 

for these objects. The light solid line in the upper left denotes the T c oc p c relation 
expected during the Kelvin-Helmholtz contraction for a fixed mass non-degenerate star. 
Deviations from this relation occur when electron degeneracy occurs, which is shown by 
the grey line at rj e^/fe^T = 4, roughly where the electron degeneracy has increased the 
electron pressure to twice that of an ideal electron gas. We extend the mass range down to 
M = 0.01M Q to reveal the distinction between main sequence stars and brown dwarfs. That 
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distinction becomes clearer in Figure 17 which shows the L evolution for a range of stars with 
M < 0.3M Q . Only the M > O.O8M stars asymptote at late times to a constant L, whereas 
the others continue to fade. We also exhibit the expected scaling for a contracting fully 
convective star with a constant T e s, L oc t~ 2 ^ 3 . At the D mass fraction of this calculation, 
Xd = 3 x 10 -5 , the onset of D burning provides some luminosity for a finite time, causing 
the evident kink at early times. 
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Fig. 16. — Trajectories of central conditions for fully convective M < O.3M stars as they 
approach the main sequence (M > O.O8M ) or become brown dwarfs for Y = 0.275 and 
Z = 0.019. Each solid line shows T c and p c for a fixed mass, M, noted at the end of the line 
(when the age is 3 Gyr). The dashed blue lines are isochrones for ages of 10, 30,100, 300 
Myr and 1 Gyr. The purple squares and red circles show where D and 7 Li is depleted by a 
factor of 100. The green triangles show the ZAMS. 
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Fig. 17. — Luminosity evolution for fully convective M < O.3M stars as they approach the 
main sequence (M > O.O8M ) or become brown dwarfs for Y = 0.275 and Z = 0.019. From 
top to bottom, the lines are for M = 0.3, 0.2, 0.15, 0.1, 0.08, 0.07, 0.05, 0.04, 0.03, 0.02, 
0.015, 0.010, 0.005, 0.003, 0.002, & 0.001M Q . The purple squares (red circles) denote where 
D ( 7 Li) is depleted by a factor of 100. 



MESA star models evolved from the PMS Hayashi line to an age of 10 Gyr with masses 
ranging from 0.09 to O.OO1M are compared with the models of Baraffe et al. (2003 , BCBAH) 
in Figure 18 Ages in increasing powers of 10 are marked by filled circles along each track 
from 1 Myr to 10 Gyr. For comparison, separate points from the BCBAH evolutionary 
models are plotted as plus symbols ("+")• So as to match the choice of BCBAH, we set 
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the D mass fraction at = 2 x 10 5 (Chabrier et al. 2000). Evolution at the youngest 



ages is uncertain due to different assumptions regarding D burning but beyond 10 Myr the 
MESA star and the BCBAH models overlap at almost every point. Note that the BCBAH 



models were only evolved to 5 Gyr for the two lowest masses shown in Figure 18 
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Fig. 18. — Evolution of very low mass stars and substellar objects from 0.09 to 0.001 M for 
Z = 0.02, Y = 0.28 in the log(g)-T e fj plane. The solid lines are the MESA star tracks, with 
labeled masses (in M ) at the bottom of each. The filled circles denote the location of each 
track at a given age. The plus symbols (+) mark the locations of the BCBAH tracks for the 
same masses and ages. The two lowest mass tracks from BCBAH do not extend to 10 Gyr. 
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7.1.2. Code comparisons of0.8M G and 1M models 



The Stellar Code Calibration Project (Weiss et al. 2007) was created to provide in- 



sight into the consistency of results obtained from different state-of-the-art stellar evo- 
lution codes. The contributors performed a series of stellar evolution calculations with 
the physics choices held constant to the greatest extent possible. This section compares 
MESA star models with published results from that project for two specific cases. The com- 



parison codes are BaSTI /FRANEC ( |Pietrinferni et al.||2004| ), DSEP ( potter et aL||2007ft , 
and GARSTEC QWeiss k Schlattl||2008[ ). MESA star models lie within the range exhibited 
by BaSTI /FRANEC, DSEP, and GARSTEC in these comparisons. 

Two examples are shown here, a O.8M , Z = 10~ 4 star and a 1M & ,Z = 0.02, both 
modeled from the pre-MS to the onset of the He core flash. The models assume, as much 



as possible, the same nuclear reaction rates (NACRE), opacities (OPAL and Alexander & 
Ferguson 1994), equation of state (FreeEOS), and mixing length (omlt = 1-6). These tests 
do not represent the best models for the various codes. Instead, the goal of the comparisons 
was to see how consistent the codes would be when using simple assumptions and comparable 
input physics ( Weiss et al.||2007 ). While the agreement is good in most respects, in temporal 
resolution there is a discrepancy. 



M=0.8 Mq, Y d = 0,25, Bo = 0.0001 y 1 ' 

■ MESA / 

DSEP / 
-- GARSTEC / 
" --■ BASTI / 










M=1.0 M a> Y = CI.28, Z = 0.02 



Io g (g cm 3 ) 




log Pc (g cm 3 ) 



Fig. 19.— Stellar Code Calibration project models for the O.8M ,Z = 10 -4 (left) and 
the 1M Q ,Z = 0.02 (right) cases. The upper-left panels show the H-R diagram; upper- 
right panels show luminosity versus central temperature; lower-left panels show central T-p; 
lower-right shows luminosity versus age. 



The H-R diagram of the 0.8M & ,Z = 10 4 case is nearly identical for all four tracks 
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except near the main sequence turnoff, where DSEP and BaSTI/FRANEC are hotter than 
MESA star and GARSTEC. These models have essentially no convection during the main 
sequence and there is remarkably little scatter during this phase. In the T c -L plane, the 
models are almost indistinguishable until they enter the red giant phase at L pa 1OL , where 
the central temperatures differ slightly only to re-converge at maximum luminosity on the 
red giant branch. Finally, the lifetime-luminosity plane indicates that the four codes split 
into two pairs with one pair shorter lived by pa 5% than the other pair. It is beyond our 
scope here to explain the reasons for these differences; the purpose of the present comparison 
is to indicate that MESA star produces results that are consistent with the range exhibited 
among the other three codes. 

Convection plays a more prominent role in the 1M , Z = 0.02 case and the scatter is 
greater than in the Z = 10 -4 case. The BaSTI/FRANEC model is hotter than the other 
three models. Treatment of convection and, in particular, the resolution of the surface 
convection zone is primarily responsible for the spread seen in the main sequence portion of 
the tracks. 

In both cases, the central conditions are very similar until the models become red 
giants. In the Z = 0.02 case, the range of lifetimes is somewhat reduced compared to the 
Z = 10~ 4 case with BaSTI/FRANEC and MESA star shortest (though the order is reversed 
with respect to the Z = 10~ 4 case) followed by DSEP and then GARSTEC. 



7.1.3. The MESA star Solar model 



MESA star performs a Solar model calibration by iterating on the difference between 



the final model and the adopted Solar parameters of L & and R Q (from Bahcall et al. 2005 ) 



and the surface value of Z s /X s from Grevesse & Sauval ( |1998 ) at 4.57 Gyr. This is done 
by iteratively varying «mlt and the initial and Z± values (all for the abundance ratios of 



Grevesse & Sauval 1998), while including diffusion. 



The properties of the converged model (which reaches the desired parameters to better 
than one part in 10 5 ) are shown in Table 10, and match the measured depth, Rqz, of the 



surface convection zone within 1-a and the surface Helium abundance, Y s , within 2-a (Bahcall 



et al. 2005). The difference between the model and the helioseismologically inferred Solar 



sound speed profile is compared with similar results from Bahcall et al. (1998, BBP98) and 



Serenelli et al. (2009, S09) in Figure 20 demonstrating that MESA star is capable of stellar 
evolution calculations at the level of 1 part in 10 3 demonstrated by others (Basu & Antia 
20081. 
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Fig. 20. — Comparisons of the sound speed profiles within the sun. The red solid line shows 
the relative difference in the sound speed between MESA star predictions and the inferred 
sound speed profile from helioseismic data (taken from Bahcall et al. (1998)). The green- 
dashed and blue-dotted lines show the same for the standard Solar models of IBahcall et al.l 
(p98| BBP98) andlSerenelli et al.| p009} S09), respectively. 



7.2. Intermediate Mass Structure and Evolution 



MESA star can calculate the evolution of intermediate mass stars ( 2 < M/M < 10) 
through the He-core burning phase and the advanced He-shell burning Asymptotic Giant 
Branch (AGB) phase. MESA star produces results compatible with published results from 
existing stellar evolution codes. 
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Table 10. MESA star Standard Solar Model at 4.57 Gyr 



Quantity Value 



Converged Input Parameters 

«mlt 1.9179113764 
Yi 0.2744267987 
Zi 0.0191292323 

Properties of Converged Model 

(Z/X) s 0.02293 

X s 0.73973 

Y s 0.24331 

Z s 0.01696 



X c 0.33550 



Z c 0.02125 

Rcz/Rq 0.71398 

logp c 2.18644 

logP c 17.3695 

logT c 7.19518 

RMS[(c M odci -c )/c ] 0.00093 
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Fig. 21. — Top: MESA star H-R diagram for 2-10 M Q models from the pre-main sequence 
to the end of the thermally pulsating AGB. Bottom: trajectories of the central conditions. 
The filled red points show the ZAMS. 



We start by showing in Figure 21 a grid of MESA star evolutionary tracks with masses 
ranging from 2 to 1OM with Z = 0.02. The top panel shows the evolution in the H-R 
diagram while the bottom panel shows the evolution in the T c — p c diagram. The 8 and 
10M Q models start to ignite carbon burning off center, whereas the 2 — 7M Q models produce 
C/O white dwarfs. The lack of a complete treatment in MESA star of liquid diffusion inhibits 
our ability to verify the resulting white dwarf cooling sequences from MESA star at this time. 
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7.2.1. Comparison of EVOL and MESA star 



We compare = 2M , Z = 0.01 stellar models from MESA star and EVOL (Blocker 



1995 Herwig 2004 Herwig & Austin 2004) starting from the pre-main sequence to the 



tip of the thermal pulse AGB (TP-AGB). Both codes employed the exponentially-decaying 



overshoot mixing treatment described by Herwig (2000, see £5.2) at all convective boundaries 
with / = 0.014, except during the third dredge-up where we adopt / = 0.126 at the bottom 
of the convective envelope to account for the formation of a 13 C pocket, and at the bottom 
of the He-shell flash convection zone we use / 



0.008 (Herwig 2005). 



In both codes we use the mass loss formula of Blocker (1995, see £6.6). Thermal pulses 



start at a slightly lower core mass, and hence luminosity, in the EVOL model. In order to 
maintain similar envelope mass evolution through the TP-AGB, the parameter t)b\ in the 
mass loss formula was set to 0.05 in MESA star and 0.1 in EVOL. Every effort has been 
made to tailor the MESA star model to the EVOL model. However, the AGB evolution is 
very sensitive to the initial core mass, which depends on the mixing assumptions and their 
numerical implementation during the preceding He-core burning phase. Consequently, small 
differences on the TP-AGB are unavoidable when comparing tracks from two codes. 



As shown in Figure 22, the EVOL and MESA star tracks compare well in the H-R 



diagram. Table 11 shows that key properties differ by less than 5%. MESA star has the 
ability to impose a minimum size on convection zones below which overshoot mixing is 
ignored. EVOL does not have such limits, leading to more mixing of He into the core and, 
hence, the ~ 4% larger age of the EVOL sequence at the first thermal pulse. 

The thermal-pulse AGB (TP-AGB) is characterized by recurrent thermonuclear insta- 
bilities of the He-shell, leading to complex mixing and nucleosynthesis. These processes are 
properly represented in MESA star calculations, as revealed in Figure 23 The ability of 



MESA star to calculate the evolution of stellar parameters in a smooth and continuous man- 
ner even during the advanced thermal pulse phases and beyond is demonstrated in Figure 
24 The top panel shows the evolution in the H-R diagram, whereas the bottom panel shows 
the evolution of the conditions in the C/O core. The adiabatic cooling in the C/O core that 
occurs during the He flash (due to the pressure dropping at the surface of the C/O core) is 
evident in the downturns that are parallel to the line of constant degeneracy (which is also 
the adiabatic slope). The overall trend of increasing p c reflects the growing C/O core mass, 



which for this model is shown in the top panel of Figure 23 
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Table 11. Comparison of MESA star and EVOL models with M t = 2M , Z = 0.01 



MESA star EVOL 



Main sequence lifetime (Gyr) 


0.939 


0.962 


Deepest penetration of first dredge-up (M ) 


0.328 


0.327 


H-free core mass at the end of He-core burning (M ) 


0.466 


0.454 


Core mass at first thermal pulse (M Q ) 


0.504 


0.481 


Age at first thermal pulse (Gyr) 


1.269 


1.328 


Core mass at 2nd thermal pulse with DUP (M ) 


0.563 


0.563 


following interpulse time (1000 yr) 


116 


106 


following pulse-to-pulse core growth (1O~ 3 M ) 


6.4 


6.9 


dredge-up mass at following pulse (1O~ 3 M ) 


1.1 


1.3 




Fig. 22. — The 2M , Z = 0.01 tracks up to the first thermal pulse from EVOL (solid black 
line) and MESA star (thick grey line) in the H-R diagram. 
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Fig. 23. — Properties of a Mi = 2M star from MESA star as it approaches the end of the 
AGB. Top: the boundaries of the C/O core and the He layer. Middle: luminosities from 
hydrogen and helium burning. Bottom: central temperature evolution. 
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log p c (g cm 3 ) 



Fig. 24. — Top: H-R diagram for the 2M MESA star model during the thermal pulses on 
the AGB. Bottom: trajectories of the central conditions in the C/O core during the thermal 
pulses. The line showing constant degeneracy is marked. 



An example of the evolution of convection zones, shell burning and total luminosities 

CIS 9b 



as well as core boundaries for two subsequent thermal pulses is shown in Figure 25 



function of model number; compare to Figure 3 in Herwig (2005 ). Quantitative comparison of 
interpulse time, core growth and dredge-up amount (see Table 11) shows excellent agreement 
between the MESA star thermal pulses and the equivalent pulses in the EVOL sequence. 
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Fig. 25. — Kippenhahn diagram with luminosities for the 2 nd and 3 rd thermal pulses with 



third dredge-up of the 2M , Z = 0.01 MESA star track shown in Figure 23 . 
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Fig. 26. — The C/O number ratio (top panel) and stellar mass, M (bottom panel) as a 
function of time from EVOL (solid black line) and MESA star (thick grey line). Time has 
been set to zero for both tracks at the onset of the third dredge-up. 



Another important property of the He-shell flashes is the intershell abundance as a result 
of the convective mixing and burning. Again, the comparison of results from both codes 
shows good agreement, which is expected since they both implement the same overshooting 
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mixing assumptions for the He-shell flash convection zone. A consequence of the third dredge- 
up is the gradual increase of the envelope C/O ratio as thermal pulses repeatedly occur. 
Since the 2M models do not experience hot-bottom burning, the evolution of this ratio is 
an effective probe of the cumulative efficiency of the third dredge-up in these simulations 



2 r , 



The top panel of Figure 26 shows the surface C/O ratio evolution according to EVOL 
(dashed-red line) and MESA star (solid black line). They are in good agreement, e.g. in 
terms of the time period over which the third dredge-up occurs and the amount by which 
C/O increases. The mass loss history over the same time period, shown in the bottom panel 
of Figure 26 , is similar by design. 



7.2.2. Interior structure of Slowly Pulsating B Stars and Beta Cepheids 



The advent of space-based asteroseismology for main sequence B stars with the Corot 
(Degroote et al. 2009[ ) and Kepler (Gilliland et al. 2010) satellites is probing the slowly 
pulsating B stars (SPBs, M « 3 — 8M ) and the more massive (M « 7 — 2OM ) (5 Cepheids 
(Degroote et al. 2009, 2010). These stars are all undergoing main sequence H burning 



and are unstably pulsating due to the k mechanism from the Fe-group opacity bump at 
2 x 10 5 K (Dziembowski et al. 1993). The observed modes have finite amplitudes 



T 



deep in the stellar core, demanding a full interior model for mode frequency (and stability) 



prediction (Dziembowski et al. 1993 Pamyatnykh et al. 2004). 



25 For massive AGB stars MESA star shows the expected hot-bottom burning behavior, including, for 
example, the avoidance of the C-star phase for a 5M Q , Z — 0.01 stellar model track despite efficient third 
dredge-up. 
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Z = 0.015 

on I — I X c = 0.2414 
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Fig. 27. — Comparison of MESA star predictions of the Brunt- Vaisala frequency, N, to two 
cases from the literature; in both cases, the MESA star model is shown as a solid line while 
the literature values are plotted as filled green circles. Comparisons are made at fixed X c for 
H burning stars. The bottom panel shows a 4M star from Dziembowski et al. (1993), and 
the top panel shows a M = 9.858M star from Pamyatnykh et al. (2004). In keeping with 
the way the numbers are presented in these papers, the vertical axes are different in the two 
panels with the bottom one in dimensionless units of N/(3GM/ R 3 ) 1 ^ 2 and the top in cycles 
per day. 
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These papers provide a few specific models that allow a direct comparison to the 
MESA star prediction of the Brunt- Vaisala frequency 



N 2 



9 



1 d In P dlnp 



dr 



T% dr 



9 



g dlnp 

c? dr 



where c 2 s = Y\Pjp is the adiabatic sound speed, and we used hydrostatic balance, dP/dr = 
—pg. Numerically, these are obtained by interpolating the sound speed at the cell boundary, 
whereas dlnp/dr is estimated by numerical differencing and then smoothed. This method 
naturally captures the extra restoring force from composition gradients, especially relevant 
in these evolving stars that leave a He rich radiative region above the retreating convective 
core during the main sequence. 



Our first comparison is to Dziembowski et al. (1993)'s M = 4M main sequence star 
with Z = 0.02 at a time when the hydrogen abundance in the convective core is X c = 0.37. 
With no overshoot from the convective core, Dziembowski et al. (1993) found log L/L & = 
2.51 and logT efr = 4.142 whereas MESA star gives log L/L Q = 2.50 and logT efr = 4.125. 



The top panel in Figure 27 compares the MESA star results (solid line) to the values (green 



circles) from Figure 3 of Dziembowski et al. (1993). The agreement is remarkable as an 



integral test of MESA star. The bottom panel of Figure [27] is a comparison to the more 
massive M = 9.858M Q main sequence star with Z = 0.015 from Figure 5 of Pamyatnykh 
et al.| Q2004D at an age (15.7 Myr) when X c = 0.2414 with log L/L Q 
4.3553. MESA star gave log L/L Q 



3.969 and logT eff = 
4.358 and an age of 16.4 Myr at the 



3.966, logT eff 

same value of X c . These comparisons highlight the readiness of MESA star for adiabatic 
asteroseismological studies of main sequence stars. 



7.3. High Mass Stellar Structure and Evolution 



To explore MESA star's results in this mass range, models of 15M Q , 2OM , and 25M 
of solar metallicity and 1000M o of zero metallicity were evolved from the Hayashi track to 
the onset of core-collapse. Nuclear reactions are treated with the 21 isotope reaction net- 



work, inspired by the 19 isotope network in Weaver et al. (1978), that is capable of efficiently 



generating accurate nuclear energy generation rates from hydrogen burning through silicon 



burning (see £4.5). This network includes linkages for PP-I, steady-state CNO cycles, a 
standard a-chain, heavy ion reactions, and aspects of photodisintegration into 54 Fe. Atmo- 



spheres are treated as a r=2/3 Eddington gray surface as described in £5.3 Mass loss for the 



solar metallicity stars uses the combined results of Glebbeek et al. (2009); Vink, de Koter &; 



Lamers (2001); Nugis & Lamers (2009); Nieuwenhuijzen & de Jager (1990), as described in 



£6.6 These massive star models are non-rotating, use no semi-convection, employ a mixing 
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length parameter of olmlt — 1-6, and adopt /=0.01 for exponential diffusive overshoot (see 



5.2) for convective regions that are either burning hydrogen or are not burning. 

Most of this section consists of comparisons to results from other stellar evolution codes. 



However, for consistency (and completeness), we show in Figure 28 the H-R diagram and cen- 
tral condition evolution of 10 — 1OOM stars from the PMS to the end of core Helium-burning. 
Though these are stars with Z = 0.02, we turned off mass loss during this calculation so that 
the plot would be easier to read and of some pedagogical use. The tendency of T c to scale 
with p c (also a constant radiation entropy) during these stages of evolution is expected 
from hydrostatic balance with only a mildly changing mean molecular weight. The rest of 
the calculations in this section included mass-loss as described above. 
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Fig. 28. — Top: H-R diagram for 10 — 1OOM models from the PMS to the end of core 
Helium burning for Z = 0.02 but with zero mass loss. Bottom: trajectories of the central 
conditions in the T — p plane over this same evolutionary period. 
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7.3.1. 25 M Model Comparisons 



Figure 29 shows the T c — p c evolution in Mj = 25M solar metallicity models from 



MESA star, Kepler (private communication - Alex Heger), Hirschi et al. (2004), and FRANEC 



(Limongi & Chieffi 2006) from helium burning until iron-core collapse. The curves fall below 
the T c oc pi 1 scaling relation as the mean molecular weight increases due to the subsequent 
burning stages. The curves are also punctuated with non-monotonic behavior when nuclear 



fuels are first ignited in shells. Figure 29 shows that MESA star produces core evolutionary 
tracks consistent with other pre-supernova efforts. The bump in the MESA star curve around 
carbon burning is due to the development of central convection whereas the other codes do 



not (although see Figure 2 of Limongi et al. 2000). The development of a convective core 



during carbon burning depends on the carbon abundance left over from core helium burning 



(Limongi et al. 2000). 



The mass fraction profiles of the inner 2.5M of this Mj = 25M model are shown 



in Figure [30] at the onset of core collapse. At the time of these plots, the infall speed has 
reached ~ 1000 km s~ l just inside the iron core (at m = 1.5M ) and the electron fraction, 
Y e , has dropped below ~ 0.48. The oxygen shell lies at 1.88 < m/M < 2.5, the silicon 



shell between 1.61 < m/M < 1.88, and the iron core at m < 1.61M . Figure 31 shows T, 
p, S, the radial velocity, the infall timescale, and Y e of this inner 2.5M . Note the entropy 
decrements at the oxygen, silicon and iron core boundaries. 



Figure 32 summarizes the history of the inner 7M of this Mj = 25M model as a 
function of interior mass (left y-axis). Evolution is measured by the logarithm of time (in 
years) remaining until the death of the star as a supernova (x-axis), which reveals the late 
burning stages. Levels of red and blue shading indicate the magnitude of the net energy 
generation (nuclear energy generation minus neutrino losses), with red reflecting positive 
values and blue indicating negative ones. The vertical lines indicate regions that are fully 
convective. Note the appearance of a convective envelope characteristic of a red supergiant 
late during helium burning. Abundance profiles of key isotopes during the major burning 
stages are shown (right y-axis). The hydrogen core shrinks towards the end of hydrogen 
burning, and the helium core grows as helium is depleted. The total mass shrinks to about 
M = 12M due to mass loss. 
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Fig. 29. — Evolution of the central temperature and central density in solar metallicity 
Mi = 25M models from different stellar evolution codes. The locations of core helium, 
carbon, neon, oxygen, and silicon burning are labeled, as is the relation T c oc pl^ 3 . 




Fig. 30. — Mass fraction profiles of the inner 2.5M of the solar metallicity M$ = 25M 
model at the onset of core collapse. The reaction network includes links between 54 Fe, 56 Cr, 
neutrons, and protons to model aspects of photodisintegration and neutronization. 
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Fig. 31. — Profiles of T (top left), p (middle left), dimensionless entropy (bottom left), 
material speed (top right), infall timescale (middle right), and electron fraction Y e = Z/A 
(bottom right) over the inner 2.5M of the Mj = 25M star at the end of the pre-supernova 
evolution. 
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log(years before end of run) 

Fig. 32. — Kippenhafm diagram showing the full time evolution of the inner 7 M & of the 
Mi = 25M evolutionary sequence from the main sequence to the onset of core collapse. Mass 
coordinate and abundance mass fraction are labeled on the left and right y-axes, respectively. 
The shaded bar on the right indicates the net energy generation: red for positive values and 
blue for negative values. The vertical lines indicate convection. 



7.3.2. Comparison of 15, 20, and 25M Models 



Now that we have shown that the Mi = 25M MESA star models compare well to 
previous efforts at the qualitative level, we will make more detailed comparisons to other 
available results. Table 12 compares the core burning lifetimes of solar metallicity stars 
with Mi = 15,20 and 25M , from MESA star, |Hirschi et a"L] ( |2004[ ), |Woosley et~aL] ( |2002[ ), 
and Limongi et al. (2000). We define a core burning lifetime to begin when the central mass 
fraction of fuel has dropped by 0.003 from its maximum value (or onset of central convection) 
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and to terminate when the central mass fraction has dropped below 10 -4 (or the end of central 
convection). Different authors adopt different lifetime definitions, which likely contribute to 
some of the scatter. The hydrogen burning lifetimes for the 15M Q , 2OM , and 25M models 



from the different authors are within 10% percent of each other, with the Limongi et al. 



(2000) models generally having the shortest lifetimes and the Woosley et al. (2002) models 



having the longest lifetimes. There is more spread in the helium burning lifetimes, with 



MESA star models showing shorter lifetimes and Woosley et al. (2002) models having the 



longest lifetimes. The carbon burning lifetimes show agreement within 20% for the 15M 
model, but differ by factors of ~3 for the 2OM and 25M models. The neon, oxygen, and 
silicon burning lifetimes show agreement within 20% between some models, but factor of 
~5 differences in others. It is beyond the scope of this paper to put the different lifetime 
definitions on the same footing, and explore the reasons for these differences. Nevertheless, 



Table 12 suggests MESA star produces lifetimes consistent with the range of lifetimes from 
other works. 



Table 13 compares pre-supernova core masses of solar metallicity stars with Mj = 15, 20 



and 25M models from MESA star, Hirschi et al. (2004), Rauscher et al. (2002), Heger et al. 



(2000), and Limongi et al. (2000). MESA star core masses are defined as the mass interior to 



the location where the element mass fraction is 0.5. The definitions used by various authors 
may differ, contributing to scatter in the results. However, most of the scatter is probably 
due to the different mass loss prescriptions used by different authors, resulting in different 



total masses. The helium yields differ by about 10%, with the Heger et al. (2000) models 



producing less helium. There is more diversity in the C+O+Ne bulk yields, up to a factor 



of 2 for the 25M model, with the Rauscher et al. (2002) models producing the most and 



the Heger et al. (2000) models producing the least. Strikingly, the Fe core masses show 



less variations, with the Hirschi et al. (2004) models producing the heaviest cores. Table 13 



suggests MESA star produces bulk yields compatible with previous efforts. 



7.3.3. 1000M o metal-free star capabilities 

We close this section with a demonstration of MESA star's capabilities by describing the 
unlikely scenario of a purely metal-free stellar evolution of a Mj = 1000M o star. The T c — p c 
trajectory for a 1000M o , zero metallicity, zero mass loss model is shown in the left panel 
of Figure 33 The starting time point is in the lower left corner and the final model, at the 
onset of core-collapse, is in the upper right at very high values of T c and p c . Fluid elements 
in the region to the left of the red-dashed line have Ti < 4/3. When the center enters this 
region, the central portions of the star become dynamically unstable and begin to contract. 
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Table 12. Massive Star Core Burning Lifetime Comparison 



Core Burning Lifetime (years) 

Element HMM WHW LSC MESA 



Mi = 15M C 



H 


1.13 




1.11 


1.07 


1.14 


xlO 7 


He 


1.34 




1.97 


1.4 


1.25 


xlO 6 


C 


3.92 




2.03 


2.6 


4.23 


xlO 3 


Ne 


3.08 




0.732 


2.00 


3.61 







2.43 




2.58 


2.43 


4.10 




Si 


2.14 




5.01 


2.14 


0.810 


xl0~ 2 






Mi 


= 2OM 








H 


7.95 




8.13 


7.48 


8.01 


xlO 6 


He 


8.75 




11.7 


9.3 


8.10 


xlO 5 


C 


9.56 




9.76 


14.5 


13.5 


xlO 3 


Ne 


0.193 




0.599 


1.46 


0.916 







0.476 




1.25 


0.72 


0.751 




Si 


9.52 




31.5 


3.50 


3.32 


xlO" 3 






Mi 


= 25M 








H 


6.55 




6.706 


5.936 


6.38 


xlO 6 


He 


6.85 




8.395 


6.85 


6.30 


xlO 5 


C 


3.17 




5.222 


9.72 


9.07 


xlO 2 


Ne 


0.882 




0.891 


0.77 


0.202 







0.318 




0.402 


0.33 


0.402 




Si 


3.34 




2.01 


3.41 


3.10 


xl0~ 3 



References. - HMM- |Hirschi et al.| p004| ; WHW- 
Woosley et aL] ( |2002[ ); LSC- |Limongi et al. ( |2000| ); MESA-this 
paper 



Table 13. Pre-Supernovae Core Mass Comparisons 



Mass (M ) HMM RHW HLW LSC MESA 



Mi = 15M ( 



Total 


13.232 12.612 


13.55 


15 


12.81 


He 


4.168 4.163 


3.82 


4.10 


4.37 


C+O+Ne 


2.302 2.819 


1.77 


2.39 


2.27 


"Fe" 


1.514 1.452 


1.33 


1.429 


1.510 


Mi = 20 M. 


Total 


15.69 14.74 


16.31 


20 


15.50 


He 


6.21 6.13 


5.68 


5.94 


6.33 


C+O+Ne 


3.84 4.51 


2.31 


3.44 


3.77 


"Fe" 


1.75 1.46 


1.64 


1.52 


1.58 


Mi = 25M 


Total 


16.002 13.079 


18.72 


25 


15.28 


He 


8.434 8.317 


7.86 


8.01 


8.41 


C+O+Ne 


5.834 6.498 


3.11 


4.90 


5.49 


"Fe" 


1.985 1.619 


1.36 


1.527 


1.62 



References. - 


- HMM- 


Hirschi et al. 


(2004) 


Rauscher et al. 


( 


2002 


); 


HLW- 


Heger et al. 



LSC-Limongi et al. (2000): MESA-this paper 
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However, the entire star does not collapse because the infalling regions become denser and 
hotter, causing the central region to leave the I 1 ! < 4/3 region and the infall to slow. Now 
another part of the star moves into the T\ < 4/3 region and begins to infall at high velocity. 
The net result is that the region where F± < 4/3 starts at the center and moves outward. 



The right panel of Figure 33 shows the material speed and Ti profiles for the final model, 
where the infalling region is now at m ~ 48OM . 



60 

o 




9 5- 




logp c (gem 3 ) 



m(M Q ) 



Fig. 33. — Time history (left panel) of T c and p c in a 1OOOM , zero metallicity, zero mass 
loss model. Also shown are the boundaries within which Ti <4/3. Material speed and T\ 
profiles (right panel) for the final model. 



The global history of the 1000M Q model as a function of time is shown in the left panel 
of Figure 34 A convective envelope appears during late helium burning. Abundance profiles 
of key isotopes during the major burning stages are shown (right y-axis). Note the short 
carbon burning era. At late times the core photodisintegrates to 4 He instead of creating 56 Ni 
because of the lower central densities encountered in these supermassive progenitors. This 
also partially causes the large endothermic central regions of the star. 
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7.4. Stellar Evolution with Mass Transfer 

MESA star can be used to examine how a star responds to mass loss or accretion (see 
£6.6). This opens up a large variety of possible applications, including accretion onto white 
dwarfs for classical novae and thermonuclear supernovae, mass transfer in tight stellar bina- 
ries, and learning the response of a star to sudden mass loss. We show two examples where 
MESA star's results can be compared to previous work. The first is a mass-transfer scenario 
relevant to P or b < 2 hour cataclysmic variables, and the second is the response of a neutron 
star to accretion of pure He. 
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7.4-1- Mass Transfer in a Binary 



To illustrate MESA star's ability to calculate the impact of mass loss on a star, we 
model the evolution of a compact binary consisting of a Roche Lobe filling low-mass ZAMS 
(M < O.2M ) model and an accreting white dwarf with Mwd — O.6M . These short orbital 
period (P rb < 2 hr) cataclysmic variables are the end points of these mass transferring 
systems QPatterson|1998j|Kolb fc Baraffe|1999) and are now being discovered in large numbers 



in the SDSS database (more than 100 studied by Gansicke et al. 2009). 



We model the parameters of the binary system and the Roche lobe overflow triggered 
mass transfer rate M as in 



Madhusudhan et al 



(2006). So as to compare to the previous 



work of Kolb & Baraffe (1999), we presume angular momentum losses from gravitational 
wave emission and keep the accreting WD mass fixed at its initial value, Mwd = O.6M . 
The evolution of the donor star is carried out by MESA star, using the r = 100 atmosphere 



tables from atm. The evolution shown in Figure 35 is followed for over 6 Gyr until the 
donor has been reduced to a brown dwarf remnant of M ~ O.O3M (see Table 14). During 



that time, the binary period drops to a minimum value of 67.4 minutes and then increases, 
independent of the initial donor mass. This plot is very similar to Figure 1 of Kolb Baraffe 



(1999). We also show in Table 14 the evolution in time of the main properties of the donor 
star and mass transfer rate of the Mj = O.21M model. Again, this agrees with the results 



in Table 2 of Kolb & Baraffe (1999). The prime differences can be attributed to a slightly 



different R(M) relation. 
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Fig. 35. — Mass transfer rate for cataclysmic variables with low mass main sequence donor 
stars of varying initial masses Mj. Each line shows the M history for different initial mass 
donors, all accreting onto a Mwd = O.6M white dwarf. After a period of initial adjustment 
to the mass transfer, each track tends to the same trajectory, showing the orbital period 
minimum at P or b = 67.4 minutes. 
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Table 14. Mass Transfer History for M< = O.21M and M WD = 0.6M ( 



Time(Gyr) P orb (hr) M/M & T off (K) log(L/L ) i?/i? logM 
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7.4-2. Rapid Helium Accretion onto a Neutron Star 



The outer envelope of an accreting neutron star is modeled in MESA star by using non- 
zero boundary conditions M c and L c (see discussion in £6.2) at a finite radius R c . This 



allows for a time dependent calculation of the thermonuclear instability that yields Type I 



X-ray bursts (Strohmayer & Bildsten 2006) for those accretion rates where the burning is 



thermally unstable (M < 10 8 M© yr x ). Such calculations have been performed with the 
KEPLER code ( [Woosley et aL]|2004l |Cyburt et alT]|2010| ) and prove very valuable in direct 
comparisons to observed Type I X-ray burst recurrence times and light curves, especially for 



the H-rich accreting "clocked burster" GS 1826-24 (Heger et al. 2007). We focus here on 



pure He accretion, relevant to neutron stars in ultra-compact binaries, such as 4U 1820-30 



(Cumming 2003). 



1.4M , R c = 10 km, L c = 3.6 x 10 34 ergs s" 



14 cm s 2 (correcting for the gravitational redshift). The initial model 



For these simulations we set M, 
and g = 2.39 x 10 

consisted of 3 x 10 25 g of pure 56 Fe and accreted pure He at M = 3 x 1O~ 9 M yr -1 . We 
require a slightly higher value of core luminosity L c /M « 0.19 keV nucleon -1 to reach the 
same ignition column depth (5 x 10 8 g cm -2 ) as 
the nuclear reaction network, including the 12 
and elements ( 23 Na, 27 A1, 31 P, 35 C1, and 39 K) that can appear as intermediates in (a, p)(p, 7) 



Weinberg et al. (2006). We used 31 species in 
C bypass reaction chain 12 C(p,7) 13 N(a,p) 16 



reactions and serve as the proton source for the C bypass (Weinberg et al. 2006) 



Figure [36] shows a snapshot of the time history of the helium burning luminosity, Ln e , 
which is periodic at the Type I burst recurrence time of 9.56 hours. This luminosity, as well 



as L, very quickly exceeds L^dd, i n which case we allow for mass loss via a wind (Paczynski 



& Proszynski 1986). We arbitrarily set our time coordinate to zero at the time of maximum 



luminosity, L, in the second burst after the start of accretion. The peak for Ln e is at 
t = -0.0269 s, and L > L Edd for the time interval -0.0047 < t < 1.2169 seconds.. 



2(i 



Figure 37 shows the evolving temperature profile during the convective burning runaway, 
where time increases upwards. Though not for the same ignition depth, this plot is very 



similar to Figure 2 of Weinberg et al. (2006), including the evolution of the location of the 



top of the convective zone (open squares). Weinberg et al. (2006) discussed in detail the 



onset of heat transport in the outer, thin, radiative layer that allows for the retreat of the 
top of the convective zone. This MESA star result is the first numerical confirmation of this 
transition for a pure helium accretor and demonstrates our ability to obtain excellent time 



2(i 



A movie of this flash (made with PGstar, see ^6.9) is at 



http : //mesa. sourcef orge . net/pdf s/nshe . 



mov 
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and mass resolution as shown in Figure [38} By using nonzero center boundary conditions so 



that the dq variables (see £6.6) cover only the relatively small envelope mass, we reach a mass 
resolution of ~ 1.5 x 1O _2O M . The timestep adjustment algorithms (£6.4) provide a smooth 
change from timesteps of almost an hour between bursts down to millisecond steps at peak 



luminosity (see middle panel in Figure 38). The secular increase in the number of cells is to 
track the accumulation of the pile of ashes from each burst. The evolution of abundances 



at the base of the convective zone is shown in Figure [39] and exhibits the presence of the 
isotopes 35 C1 and 39 K. 
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Fig. 36. — The helium-burning luminosity, Lu e , as a function of time for a neutron star of 
mass M c = 1.4M and radius R c = 10 km accreting pure helium at M = 3 x 1O _9 M yr _1 . 
The Type I X-ray bursts occur every 9.56 hours. 
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log column depth (g cm 2 ) 



Fig. 37. — The evolving temperature profile during the convective burning phase of a Type 
I burst, as a function of column depth, P/g. Starting from the bottom, each successive solid 
line is the temperature profile at a later time. The open squares marks the location of the 
top of the convective zone. The top curve is at t — —0.00716. 
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Fig. 38. — The He burning luminosity, timestep, and number of cells as a function of model 
number for the MESA star simulation of an accreting neutron star. The timestep ranges 
from a millisecond to an hour, whereas the number of cells only grows by ~ 25% during the 
burst and shows a secular trend upward as partially burned material accumulates. 
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Fig. 39. — Abundances of the dominant isotopes at the base of the convective zone as a 
function of time during the Type I burst. The temperature at the base of the convective 
zone at t = 1 second is logT = 9.15. 



We have performed simulations at lower accretion rates but these become dynami- 
cal events where the temperature rises on a local dynamical timescale and are beyond the 



present scope of MESA star. While multi-dimensional hydrodynamical codes (e.g., Zingale 
eTaLl[2009| may be needed to follow the details of such an event, MESA star can be used for 
studying the longer timescale, hydrostatic evolution leading up to the point where hydrody- 
namic effects become dominant. 
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8. Summary and conclusion 

Modules for Experiments in Stellar Astrophysics (MESA) provides open source, portable, 
robust, efficient, thread-safe libraries for stellar astrophysics and stellar evolution. It provides 
tools for a broad community of astrophysicists to explore a wide range of stellar masses and 
metallicities. State-of-the-art modules include the equation of state, opacity, nuclear reaction 
rates and networks, atmosphere boundary conditions, and element diffusion. MESA features 
a modern code architecture and run-time environment. 

MESA star solves the fully coupled structure and composition equations simultaneously 
and is capable of calculating full evolutionary tracks without user intervention. It implements 
adaptive mesh refinement, sophisticated timestep adjustment, mass loss and accretion, and 
parallelism based on OpenMP. 

MESA is subjected to an ongoing testing and verification process. Current capabilities in- 
clude evolutionary tracks of very low mass stellar objects and gas giant planets, intermediate 
mass stars, pulsating stars, accreting compact objects, and massive stars from the pre-main 
sequence to late times. Future versions of MESA will include the addition of a variety of new 
physics modules, features driven by the MESA user community, and architectural refinements. 
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erer, Marco Limongi, Marcin Mackiewicz, Georgios Magkotsios, Lars Mattsson, Dan Meiron, 
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A. Manifesto 

MESA was developed through the concerted efforts of the lead author over a six year 
period with the engagement and deep involvement of many theoretical and computational 
astrophysicists. The public availability of MESA will serve education, scientific research, and 
outreach. This appendix describes the scientific motivation for MESA, the philosophy and rules 
of use for MESA, and the path forward on stewardship of MESA and advanced development 
of future research and education tools. We make MESA openly available with the hope that 
it will grow into a community resource. We therefore consider it important to explain the 
guiding principles for using and contributing to MESA. Our goal is to assure the greatest 
usefulness for the largest number of research and educational projects. 

A.l. Motivation for a new tool 

Stellar evolution calculations (i.e. stellar evolution tracks and detailed information about 
the evolution of internal and global properties) are a basic tool that enable a broad range 
of research in astrophysics. Areas that critically depend on high-fidelity and modern stel- 
lar evolution include asteroseismology, nuclear astrophysics, galactic chemical evolution and 
population synthesis, compact objects, supernovae, stellar populations, stellar hydrodynam- 
ics, and stellar activity. New observational capabilities are emerging in these fields that place 
a high demand on exploration of stellar dependencies on metallicity and age. So, even though 
one dimensional stellar evolution is a mature discipline, we continue to ask new questions 
of stars. The emergence of demand requires the construction of a general, modern stellar 
evolution code that combines the following advantages: 

• Openess: should be open to any researcher, both to advance the pace of scientific 
discovery, but also to share the load of updating physics, fine-tuning, and further 
development. 

• Modularity: should provide independent, reusable modules. 

• Wide Applicability: should be capable of calculating the evolution of stars in a 
wide range of environments, including low- and massive stars, binaries, accreting, mass- 
losing stars, early and advanced phases of evolution etc. This will enable multi-problem, 
multi-object physics validation. 
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• Modern techniques: should employ modern numerical approaches, including high- 
order interpolation schemes, advanced AMR, simultaneous operator solution; should 
support well defined interfaces for related applications, e.g. atmospheres, wind simu- 
lations, nucleosynthesis simulations, and hydrodynamics. 

• Microphysics: should allow for up-to-date, wide-ranging, flexible and modular micro- 
physics. 

• Performance: should parallelize on present and future shared-memory, multi-core/thread 
and possibly hybrid architectures so that performance continues to grow within the new 
computational paradigm. 

A tool that combines the above features is a significant research and education resource 
for stellar astrophysics. We acknowledge that some important aspects of stars are truly three- 
dimensional, such as convection, rotation, and magnetism. Those applications remain in the 
realm of research frontiers with evolving understanding and insights, quite often profound. 
However, much remains to be gained scientifically (and pedagogically) by accurate one- 
dimensional calculations, and this is the present focus of MESA. 

A. 2. MESA philosophy 

The MESA code library project is open. It explicitly invites participation from anybody 
(researchers, students, interested amateurs). Participation in MESA can take a wide range of 
forms, from just using a MESA release for a science project, to testing and debugging (i.e. 
report bugs, find fixes and submit them for inclusion into the next release) as well as taking 
on responsibility for the continued stewardship of certain aspects (modules) of the code. The 
participation of experienced stellar evolution experts is very welcome. 

Users are encouraged to add to the capabilities of MESA, which will remain a community 
resource. However, use of MESA requires adherence to the "MESA code of conduct" : 

• That all publications and presentations (research, educational, or outreach) deriving 
from the use of MESA acknowledge the Paxton et al. (2010) publication and MESA 
website. 

• That user modifications and additions are given back to the community. 

• That users alert the MESA Council (see below) about their publications, either pre- 
release or at the time of publication. 
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• That users make available in a timely fashion (e.g., online at the MESA website) all 
information needed for others to recreate their MESA results - "open know how" to 
match "open source." 

• That users agree to help others learn MESA, giving back as the project progresses. 
Users are requested to identify themselves by name, email contact, and home location. 

A. 3. Establishment of the MESA council 

The MESA project began as an initiative to construct a reliable computational tool for 
stellar structure and evolution that takes full advantage of modern processor architectures, 
algorithms and community engagement. The release of MESA has forced some explicit thinking 
of what structure is needed so as to achieve the mission of stewarding MESA in its use for 
scientific research, education and outreach, while also enabling the development of new tools 
and ideas. The MESA operating principles are simple: be open in your scientific discussions, 
give credit to all contributors, and be prepared to give back to the community of users. We 
hope that this creates an environment where the young are encouraged to become engaged 
in a career-enhancing manner. 

We have established the MESA Council that consists of those engaged in working towards 
the shared missions outlined here: 

• Steward MESA There are many ways this will be done: supporting the contributors, 
maintaining the web access and web page updates, seeking enabling funding, hold- 
ing yearly working groups that allow for continued engagement, documenting MESA 
development in the refereed literature, and sustaining advanced development. 

• Interface with the User Community This starts with answering questions from 
users, developing a way to accept new code in an integrated fashion, maintain a user 
registry, and identify new MESA Council members from those most active and engaged 
in the intelligent use of MESA. 

• Enable Scientific Research and Education with MESA Promote MESA and its 
goals, e.g., through scientific contributions at relevant conferences. Identify science 
opportunities that match MESA capabilities and facilitate and encourage appropriate 
collaborative activities. Track the science carried out by the community with MESA. 
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B. Code testing and verification 

An important part of the ongoing development of a large, complex software project, such 
as MESA, is regular, systematic testing. Testing is necessary to ensure that MESA continues 
to function as expected and that the addition of new features does not have unintended 
consequences for existing features. 

MESA is tested at the module level each time it is compiled from the install scripts. These 
tests check that each module produces results that are consistent with expectations. The 
next level is the MESA star test suite, which consists of various evolutionary cases that are 
intended to cover a broad range of applications, including Roche lobe overflow, the He core 
flash in a low mass star, the evolution of sub-stellar mass objects, advanced nuclear burning 
in massive stars, accreting white dwarfs and neutron star envelopes, and more are being 
added all the time. The test cases come in both short and long varieties. Run in serial, 
the full set of short tests completes in less than one hour on modern hardware. The long 
tests might each take one or several hours to complete. Many of the evolutionary sequences 
presented in £j7] are included in the test suite. Short tests include the very low mass models 



evolved to 10 Gyr (Figure 18) and the O.8M , Z = 10 4 track (Figure 19) while longer tests 



include the Solar model calibration (Figure 20), the "hands off" 1M Q pre-main sequence to 



white dwarf calculation (Figure 13), and Si-burning in a 15M Q} Z = 0.02 model (Tables 12 



and 13). 



The test suite is readily extended in order to ensure regular testing of certain aspects 
of MESA that are not covered by the existing set but are important for a particular avenue of 
research. A template is provided to encourage the creation of new test cases. 
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